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ABSTRACT 


Principles  of  discrete  Fourier  and  spectral  analysis 
have  been  applied  to  the  motion  of  optically  detectable 
particles  in  nerve  axons.  Particle  positions  recorded  by 
means  of  cinemicrophotography  were  digitized  from 
consecutive  series  of  film  frames.  The  trends  in  x  and  y 

position  with  time  ( whose _ values  were  the  particle  mean 

velocity  components  Vx  and  V y  respectively)  were  removed 
from  the  data  and  the  remaining  deviations  Fourier  analyzed. 
Power  spectra  obtained  from  these  analyses  indicated  that 
the  vast  majority  of  the  signal  power  in  the  deviations  fell 
below  .1  Hz.,  indicating  that  the  filming  rate  (3  fps)  had 
been  sufficiently  rapid  to  record  all  but  the  most  minute 
details  of  particle  motion.  All  power  spectra  contained  a 
strong  peak  at  low  frequencies,  and  it  was  concluded  that 
variations  in  the  instantaneous  velocities  of  the  particles 
could  be  character ised  as  low  frequency  oscillations. 

The  trend -free  measurements  of  particle  x  and  y 
positions  were  processed  with  a  band- limited  dif ferent iator 
to  obtain  estimates  of  the  instantaneous  velocity 
components.  The  character ist ics  of  the  di f ferent iator  and 
the  means  of  implementing  it  were  based  on  studies  of 
simulated  particle  motion  and  on  the  power  spectra  of  the 
motion  of  real  particles. 

Some  of  the  errors  inherent  in  the  velocity  estimation 
procedure  were  assessed  by  redigitizing  the  trajectory  of  a 
single  particle  several  times.  This  approach  yielded 
estimates  of  error  in  the  velocity  estimates  which  had  a 
standard  deviation  of  about  .02  um/sec,  or  approximately  one 
order  of  magnitude  less  than  that  of  the  finite  difference 
1  method  of  estimating  instantaneous  velocity  (as  used  by 
other  workers). 

One  hundred  seventy  two  particles  tracked  in  ten  axons 
were  analyzed.  The  x  component  of  velocity  (longitudinal  to 
the  axon)  was  typically  one  to  two  orders  of  magnitude 
greater  than  that  of  the  y  component  (radial  to  the  axon). 
The  estimates  of  the  x  component  of  instantaneous  velocity 
were  quantitatively  examined,  and  in ter -group  comparisons 
made  of  a  number  of  parameters  of  particle  motion.  The 
spatial  characteristics  of  the  velocity-position  data  were 
investigated  with  the  use  of  the  SURFACE  II  spatial  analysis 
and  display  system.  Contour  plots  of  the  x  component  of 
particle  instantaneous  velocity  as  a  function  of  position 
were  produced.  Some  correlation  between  particle  locations 
and  corresponding  instantaneous  velocities  was  apparent. 
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CHAPTER  1 
INTRODUCTION 


Each  neuron  has  a  long  cytoplasmic  extension  known  as 
an  axon  (figure  1.1).  The  axon  and  nerve  endings  have 
essentially  no  means  for  protein  synthesis,  hence  proteins 
necessary  for  the  viability  of  the  axon  must  be  obtained 
from  outside  the  axon  (Lasek, 1970) .  It  has  been  shown  that 
the  neuron  cell  body  synthesizes  most  of  the  axonal 
proteins,  and  that  these  are  moved  along  the  axon  by  a 
process  known  as  axoplasmic  transport  (Lasek, 1970; 
Ochs, 1972;  Grafstein , 1975) . 

Normal  axoplasmic  transport  is  essential  for 
maintenance  of  the  axon;  unusual  or  disturbed  patterns  of 
axoplasmic  transport  have  been  associated  with  a  number  of 
neurological  disorders  (Ochs, 1975).  If  the  motive  mechanisms 
underlying  normal  axoplasmic  transport  could  be  understood, 
it  is  probable  that  this  knowledge  would  contribute  to 
treatment  and  possible  prevention  of  these  disorders. 
Researchers  have  proposed  a  number  of  models  for  the  motive 
mechanism  but  to  date  little  evidence  supporting  any  of 
these  has  been  presented  ( Schmi tt , 1969;  Ochs, 1971; 
Lasek, 1970).  Clearly,  detailed  quantitative  analysis  of  the 
axoplasmic  transport  process  is  necessary  to  provide  a  basis 
for  a  model  adequately  reflecting  the  dynamics  of  particle 
mot  ion. 

To  describe  the  direction  of  movement  within  the  axon, 
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Figure  1.1  Nerve  cell  with  some  of  its  synaptic  and 
neuromuscular  contacts.  Inset  depicts  microscopic  view 
of  axon.  Adapted  from  Figure  1  of  Nerve,  Muscle ,  and  Synapse 
(Katz, 1966) . 
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two  terms  shall  be  used:  "somatopetal " ,  which  refers  to 
motion  towards  the  soma  (cell  body);  and  "somatofugal " , 
which  refers  to  motion  away  from  the  soma.  Exper indent ers 
have  shown,  using  a  variety  of  techniques,  that  the 
transport  process  consists  of  at  least  two  components:  (1)  a 
slow  somatofugal  movement  of  possibly  the  entire  axoplasmic 
column  at  approximately  1  mm/day  (Weiss  and  Hiscoe, 1948) , 
and  (2)  a  rapid  movement  of  materials  within  the  axoplasmic 
column  at  rates  of  the  order  of  100-1000  mm/day  (Lasek,1970; 
Heslop,1975;  Schwartz, 1979) .  In  addition,  these  two 
components  may  be  further  distinguished  from  one  another  by 
the  observation  that  the  slow  component  is  unidirectional 
(i.e.  from  the  cell  body  out  to,  and  along,  the  axon) 
whereas  the  fast  component  is  bidirectional  ( Bur dwood, 1965; 
Cooper  and  Smith, 1974;  Schwartz, 1979) . 

Until  quite  recently,  axoplasmic  transport  was  studied 
primarily  with  the  use  of  indirect  methods,  including  the 
rate  of  displacement  of  radioactive  isotopes  (Lasek,1968; 
Ochs  and  Johnson, 1969)  ,  and  of  cytochemical  markers 
(Lubinska, 1964;  Dahlstrom, 1967 ) .  These  methods,  however, 
have  certain  common  limitations.  For  example,  although  net 
movement  over  a  time  frame  of  hours  can  be  detected,  details 
of  the  movement  during  much  shorter  intervals  cannot  be 
resolved.  These  methods  lend  themselves  to  the  measurement 
of  movement  only  in  large  numbers  of  neurons  whose  axons  run 
together  for  an  appreciable  distance.  They  therefore  do  not 
provide  data  about  transport  in  single  axons. 


r 


4 


Light  microscopy  of  tissue-cultured  neurons  has 
revealed  numerous  intra- axonal  particles  moving  both  to  and 
from  the  cell  body  (Matsumoto, 1920;  Hughes, 1953;  Pomerat  et 
a  1., 1967).  Although  the  particles  moved  at  velocities 
comparable  to  rates  of  fast  axoplasmic  transport  observed  in 
adult  neurons  (Burdwood, 1965) ,  there  has  been  some  dispute 
as  to  whether  conclusions  concerning  axoplasmic  transport, 
reached  from  work  on  tissue-cultured  neurons,  can  be  applied 
to  normal  adult  neurons  (Weiss, 1967) .  Quite  recently 
however,  rapidly  moving  particles  with  dimensions 
approximately  that  of  the  wavelength  of  light  have  been 
detected  optically  in  living  nerve  fibres  isolated  from 
adult  vertebrates  (Smith,  1971,1972,1973;  Kirkpatrick  et 
a  1., 1972;  Kirkpatrick  and  Stern, 1973;  Cooper  and  Smith, 1974; 
Smith  and  Koles, 1976;  Hammond, 1977 ;  Forman  et  al . , 1977a, b) . 
These  workers  have  shown  that  it  is  now  possible  to  analyse 
the  movements  of  optically  detectable  particles  in  adult 
neurons,  and  to  relate  these  results  to  those  obtained  from 
t issue -cul tured  neurons,  and  to  those  results  obtained  from 
adult  nerves  using  other  techniques. 

Work  conducted  by  a  number  of  reseachers  suggests  a 
correlation  between  the  movement  of  membrane  bounded  bodies, 
among  them  the  optically  detectable  organelles,  and  the 
rapid  component  of  axoplasmic  transport  as  determined 
indirectly  with  biochemical  methods  ( Schwartz, 1979) . 
Assuming  that  such  is  the  case,  by  observing  and  analyzing 
the  movement  of  individual  particles  the  dynamics  of  fast 
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axoplasmic  transport  may  be  be  better  understood. 

Particle  trajectories  may  be  easily  represented  in  the 
Cartesian  co-ordinate  system.  The  microscope's  viewing  plane 
can  be  described  by  the  x  and  y  co-ordinates,  with  the  z  co¬ 
ordinate  in  the  plane  perpendicular  to  it.  Only  movements  in 
the  x  and  y  co-ordinates  may  be  directly  measured,  as  small 
movements  in  the  z  coordinate  are  not  measureable  by  this 
method.  Assuming  a  particle's  continuous  position  in  time; 
x(t),y(t);  is  known,  the  corresponding  velocities  Vx(t)  and 
V y(t)  can  be  calculated  from: 


Vx( t ) =  d  x(t) 
dt 


Vy ( t ) =  d  y ( t ) 
dt 


(1.1) 


The  particle's  radial  component  of  velocity  Vd(t)  relative 
to  its  initial  position  is  then  given  by: 


Since  particle  trajectories  are  commonly  recorded  on  16 
mm  film  the  position  of  a  given  particle  is  known  only  at 
discrete  points  in  time  (i.e.  each  frame).  Common  practice 
has  been  to  orient  the  longitudinal  axis  of  the  axon  such 
that  on  the  film  records  it  is  parallel  to  the  x  coordinate 
axis  of  the  position  recording  apparatus.  The  film  may  be 
viewed  frame  by  frame,  and  the  particle's  successive  x  and  y 
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positions  recorded.  The  sampled  trajectory  (in  two 
dimensions)  is  then  represented  as  a  series  of  numbers 

x(kT),y(kT)  where  k=0,1t2 .  N-1  (  N  is  the  number  of 

readings  and  T  is  the  sampling  period).  In  other  words,  the 
particle's  continuous  position  in  time;  x(t),y(t);  becomes  a 
discrete  time  series  x(kT),  y(kT),  and  is  suitable  for 
analysis. 

Although  this  approach  has  been  taken  by  a  number  of 
workers,  the  analytical  techniques  used  to  date  have  lacked 
sophistication.  It  has  been  common  practice  to  calculate  the 
mean  velocity  of  a  particle  from  the  time  series  associated 
with  its  trajectory,  according  to  the  following: 


V=  V<x((N-1)T)-x((0)T)H  +  {y((N-1)T)-v((0)T)P 

(N-1 )T 


(1.3) 


That  is,  the  mean  velocity  of  a  particle  over  a  given  period 
of  time  is  given  by  the  total  distance  moved  divided  by  the 
total  elapsed  time.  Mean  velocities  have  also  been  estimated 
from  the  slopes  of  linear  regression  lines  calculated  for 
particle  trajectory  distance-time  data.  Recently,  some 
workers  have  attempted  analysis  of  the  small  incremental 
movements  of  particles  by  calculating  the  particles' 
"instantaneous"  velocity  (Berlinrood  et  al.,1972;  Cooper  and 
Smith, 1974;  Breuer  et  al.,1975;  Forman  et  al ., 1975, 1977a, b) : 
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V((k+1)T)=  V{x( (k+1 )T)-x(kT)}2+{y( (k+1 ) T) -v(kT) } 2 

T 

T=sampl  ing  period;  k=0,1,2, .  N- 1 

(1.4) 

That  is,  the  instantaneous  velocity  was  estimated  from  the 
distance  which  a  particle  moved  between  successive  motion 
picture  frames.  The  sign  of  V((k+1)T)  was  defined  to  be 
positive  if  there  was  a  net  movement  between  samples  in  the 
positive  x  direction.  However,  experience  has  shown  that 
repeated  tracking  and  analysis  of  the  same  particle 
trajectory  using  this  method  yields  instantaneous  velocity 
estimates  which  widely  differ.  Absent  from  the  literature  is 
any  appreciation  of  such  errors  and  the  limitations  inherent 
in  the  method. 

Equation  1.4  is  a  first  order  approximation  to  an  ideal 
different iator .  Two  major  problems  are  encountered  in  the 
use  of  this  approximation  which  make  it  unsuitable  for 
estimation  of  particle  instantaneous  velocities:  (a)  this 
approximation  uses  only  the  first  two  terms  of  an  infinite 
series  of  terms  describing  an  ideal  di f ferent iator ,  and  is 
therefore  quite  inaccurate;  and  (b)  this  approximation  is 
inordinately  sensitive  to  random  components  in  the 
particles7  posi t ion- t ime  data,  and  hence  the  more  orderly 
and  interesting  components  may  be  obscured.  In  view  of  the 
foregoing,  it  is  suggested  that  the  use  of  a  more 
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appropriate  approximation  would  yield  more  accurate  and  more 
consistent  results,  and  could  therefore  shed  far  more  light 
on  the  mechanisms  of  axoplasmic  transport  than  currently 
exists.  Improvement  in  the  instantaneous  velocity  estimation 
technique  could  consist  of  the  following  procedures:  (a)  the 
use  of  spectral  estimation  techniques  to  examine  the 
frequency  components  of  particle  motion;  (b)  the  design  of  a 
more  suitable  approximation  of  an  ideal  different iator , 
based  on  the  results  of  (a);  and  (c)  optimization  of  the 
technique  through  analysis  of  simulated  and  actual 
trajector ies. 

This  thesis  outlines  the  development  of  an  improved 
technique  for  the  estimation  of  particle  instantaneous 
velocities  through  implementation  of  these  ideas.  The 
improved  technique  was  used  to  estimate  the  instantaneous 
velocity  of  a  number  of  somatofugally  and  somatopetally 
moving  particles  in  amphibian  axons.  The  results  thus 
obtained  were  examined  with  a  variety  of  methods,  including 
correlation  and  regression  analysis,  and  contour  plotting, 
in  an  effort  to  identify  patterns  in  particle  movement. 
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CHAPTER  II 


FOURIER  ANALYSIS 


2. 1  Introduct ion 

Fourier  analysis  of  a  continuous  function  of  time  is  a 
decomposition  of  the  function,  using  the  Fourier  transform, 
into  a  sum  of  sinusoidal  components.  By  using  various 
properties  of  the  Fourier  transform,  it  is  possible  to 
estimate  the  power  spectra  of  such  functions  and  perform 
operations  such  as  convolution. 

The  principles  of  Fourier  analysis  may  be  applied  in  a 
number  of  ways,  depending  on  the  nature  of  the  function 
encountered.  A  function  may  describe  a  consistent, 
non -random  mechanism  which  acts  in  an  "orderly"  manner.  It 
is  sometimes  possible  to  derive  a  model,  based  on  physical 
laws,  which  permits  exact  calculation  of  the  value  of  the 
function  at  any  instant  in  time.  If  exact  calculation  is 
possible,  such  a  model  is  said  to  be  deterministic. 
Principles  of  Fourier  analysis  may  be  used  to  filter  and 
convolve  such  functions. 

Alternatively,  a  function  may  act  in  an  unpredictable 
manner,  randomly  varying.  Nevertheless,  it  may  be  possible 
to  derive  a  model  which  can  be  used  to  calculate  the 
probability  of  a  future  value  lying  between  two  specified 
limits.  Such  a  model  is  known  as  a  stochastic  model.  The 
Fourier  transform  is  a  valuable  tool  in  the  estimation  of 
the  power  spectra  of  such  functions. 
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Most  physically -occur ring  functions  have  both 


deterministic  and  random  components,  in  varying  proportions. 
The  movement  of  particles  in  nerve  axoplasm  is  such  a  case. 
The  particles  are  moved  by  some  orderly  motive  forces,  and 
by  random  forces  which  result  in  a  superimposed  component  to 
the  motion.  Random  effects  caused  by  microscope,  camera, 
nerve,  and  projector  movement,  digitizer  noise  and  observer 
error  all  probably  contribute  to  the  apparent  random 
component  observed  in  the  movement  of  the  particles. 
Therefore,  complete  analysis  of  the  movement  of  particles  in 
nerve  axoplasm  requires  consideration  of  techniques  of 
Fourier  analysis  used  in  spectral  estimation,  and  in 
filtering  (convolution). 

2.2  The  Four ier  Transform 

The  Fourier  transform  X ( f )  of  a  continuous  function  of 
time  x(t)  is  given  by: 


x(t)e  dt 


(2.1) 


The  inverse  Fourier  t ran form  of  X(f)  is  given  by: 


X(f  )e 


df 


(2.2) 


It  is  well  known  that  the  output  c(t),  of  a  linear 
system  whose  impulse  response  is  h(t) ,  to  an  input  x(t),  is 


J- 
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given  by: 

^  oO 

c(t)=  J  x(T)h(t-T)dT 

—  oO 

(2.3) 

This  operation  between  x(t)  and  h ( t )  is  termed  convolution. 

In  the  frequency  domain,  the  Fourier  transforms  of 
c(t),  x(t),  and  h(t)  can  be  shown  to  be  related  by: 


C(f )=X(f )H(f ) 

and  the  output  in  the  time  domain  given  by: 

T  iz*ft 

c(t)=  J  X(f)H(f)e  df 

—  oo 


(2.4) 


(2.5) 


Hence,  if  a  system  frequency  response  H(f)  is  known, 
the  output  c(t)  to  any  input  x(t)  may  be  calculated  by 
Fourier  transforming  x(t),  multiplying  the  result  by  H(f), 
and  then  inverse  Fourier  transforming  the  resultant  product. 
Using  this  procedure  it  is  possible  to  perform  a  variety  of 
filtering  tasks,  including  differentiation  and  low  pass 
f i Iter ing. 

2.3  The  Discrete  Four ier  Transform  and  the  FFT  algorithm 

2.3.1  Introduct ion 

The  movement  of  optically  detectable  particles  in  axons 


. 
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Is  usually  recorded  on  film  by  means  of 
cinemicrophotography .  The  framing  rate  determines  the  number 
of  "samples"  of  particle  position  taken  in  a  given  period  of 
time  and  should  be  sufficiently  rapid  to  sample  the  major 
characteristics  of  particle  motion.  The  frame  exposure 
period  determines  the  length  of  time  over  which  each  sample 
is  taken.  It  must  be  short  enough  to  allow  the  recorded 
image  to  be  considered  representat ive  of  the  particle's 
position  at  an  instant  in  time.  Trajectories  are  recorded  by 
digitizing  the  x  and  y  coordinate  positions  of  particles 
over  series  of  consecutive  frames.  The  series  of  coordinate 
readings  thus  obtained  are  discrete  rather  than  continuous 
functions  of  time,  and  are  known  as  discrete-t ime  signals. 

In  light  of  the  foregoing,  the  Fourier  analysis 
described  in  the  previous  section  must  be  modified  to 
accomodate  the  discrete  nature  of  the  data.  The  Discrete 
Fourier  Transform  (DFT)  of  a  set  of  N  numbers  x(kT), 
k=0, 1 ,2. . .N-1 ,  is  given  by: 


N-1  -iZTT^FkT 

X ( mF )  =  1/N  2Z  x(kT)e 

k=0 


(2.6) 


and  the  inverse  transform  : 


I 

:  ;  roi  *H9  \ 
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N-1 

x(kT ) =  YZ  X(mF)e 
m=0 


j2’irrr»FkT 


(2.7) 


F=1/NT;  T=sampling  period 
for  k=0,1...N-1;  m=0,1...N-1 


If  equation  (2.6)  were  evaluated  in  the  obvious  manner, 
N2  complex  multiplications  and  additions  would  be  required. 
Although  more  efficient  algorithms  for  evaluating  the 
transform  had  been  developed  as  early  as  1942  (Danielson  and 
Lanczos, 1942) ,  it  was  not  until  the  now  classic  paper  by 
Cooley  and  Tukey  (Cooley  and  Tukey,1965)  describing  a 
computationally  very  efficient  algorithm  for  evaluating  the 
DFT,  that  the  DFT  has  found  widespread  application.  This 
algorithm  has  come  to  be  known  as  the  Fast  Fourier  Transform 
(FFT ) .  For  the  very  favourable  situation  when  N  is  equal  to 
an  integer  power  of  2,  the  FFT  algorithm  requires  only 
N  Log  a  N  operations.  Obviously,  computation  time  is 
considerably  reduced  with  the  result  that  the  DFT 
(implemented  via  the  FFT)  has  become  one  of  the  most 
important  tools  in  the  digital  computer  analysis  of  signals. 

The  FFT  algorithm  has  been  widely  applied  to:  (1)  the 
estimation  of  the  power  spectra  of  signals,  and  (2)  the 
convolution  of  two  signals,  specific  examples  of  which  are 
auto  and  crosscorrelat ion .  In  both  applications,  however, 
two  major  problems  associated  with  sampled  discrete- t ime 
signals  and  the  use  of  the  DFT  appear: 
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(1)  If  the  sampling  rate  is  not  sufficiently  high,  high 
frequency  components  in  the  continuous  signal  will 
masquerade  as  low  frequency  ones  in  the  sampled,  discrete 
time  series  representat ion.  This  phenomenon  is  known  as 
aliasing, 

(2)  In  general,  the  data  being  analyzed  is  composed  of 
arbitrary  frequencies.  However,  the  DFT  only  calculates  the 
Fourier  coefficients  at  integer  multiples  of  1/NT  Hz.  This 
gives  rise  to  the  problem  known  as  the  picket  fence  effect. 

2.3.2  Aliasing 

In  the  analysis  of  discrete  time  representations  of 
continuous  signals,  aliasing  describes  the  process  by  which 
high  frequency  components  of  a  signal  can  be  confounded  with 
low  frequencies  if  the  sampling  rate  is  not  sufficiently 
high  (fig.  2.1).  The  sampling  theorem  of  Shannon  (Shannon 
and  Weaver, 1963)  states  that  a  continuous  bandlimited  signal 
containing  power  from  0  to  P  Hz  may  be  uniquely  represented 
by  discrete  samples  spaced  at  T  second  intervals  if 
1/2T  >  P.  That  is,  the  sampling  frequency  (1/T  Hz)  must  be 
equal  to  or  greater  than  twice  the  highest  frequency 
component  in  the  original  signal.  1/2T  Hz  is  known  as  the 
folding  frequency,  since  frequencies  above  this  limit  are 
folded  (or  aliased)  into  the  frequency  interval  below  the 
folding  frequency. 

In  practice,  it  may  be  necessary  to  sample  at  a 
frequency  below  this  desired  minimum  sampling  rate.  If  the 
sampling  rate  is  chosen  such  that  the  percentage  of  signal 
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(a)  Example  of  a  sinusoid  of  frequency  Hz  masquerading  as  one  of 
lower  frequency  F|  Hz  when  the  sampling  frequency  Fs  (1/T  Hz  ) 
has  not  been  sufficiently  rapid. 


|X(f)| 
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(b)  As  Fh  is  higher  than  the  folding  frequency  (Fs/2),  the  power 


of  the  signal  appears  at  F|  Hz 


where  F.  =  Fu  -  Fc 


Figure  2.1  Aliasing 
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power  lying  above  the  folding  frequency  is  sufficiently  low, 
the  error  introduced  by  aliasing  will  be  small. 
Alternatively,  the  continuous  signal  may  be  low  pass 
filtered  (in  effect,  artificially  bandlimited)  prior  to 
samp ling. 

2.3.3  Picket  fence  effect 

The  DFT  calculates  the  Fourier  coefficients  at  discrete 
frequencies,  each  separated  from  the  other  by  1/NT  Hz.  At 
each  discrete  Fourier  frequency  there  exists  a  spectral 
window,  as  is  illustrated  in  figure  2.2(a).  (Note  that  the 
s ide lobes  are  omitted  for  clarity.)  If  the  transformed  data 
consists  of  components  at  only  the  integer  multiples  of  1/NT 
Hz,  no  error  will  occur.  If,  however,  components  exist  at 
frequencies  which  are  not  confined  to  these  discrete 
frequencies,  the  so-called  picket  fence  effect  requires 
consider at  ion. 

The  worst  error  occurs  when  a  signal  component  falls  at 
the  midpoint  between  two  coefficients,  and  the  signal  power 
as  seen  by  the  two  adjacent  spectral  windows  is  only  .406  of 
its  true  value.  The  effect  then  is  one  of  viewing  the  true 
spectrum  through  a  set  of  pickets. 

By  artificially  lengthening  the  data  record  with  zero 
terms  to  a  new  length  N'  ,  the  DFT  will  calculate  a  set  of 
Fourier  coefficients  at  integer  multiples  of  1/N'T  lying 
between  the  integer  multiples  of  1/NT  Hz.  This  procedure, 
known  as  "padding"  the  data  record,  does  not  affect  the 
width  of  the  spectral  window,  which  is  dependent  only  on  the 
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(a)  Picket  fence  effect 


(b)  Reduction  of  picket  fence  effect  through 

computation  of  additional  Fourier  coefficients. 


Figure  2.2  Picket  fence  effect 
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"true"  length  of  the  data.  The  result  is  considerably  more 
overlap  of  adjacent  windows.  This  is  demonstrated  in  figure 
2.2(b)  for  the  case  N'=2N.  The  consequences  of  the  picket 
fence  effect  are  further  reduced  in  most  cases  by  using 
non -rectangular  data  windows  with  wider  main  lobes. 

2.4  Spectral  est imat ion 

2.4.1  Introduct ion 

The  autocorrelat ion  function  R(kT)  of  an  infinite 
series  x(nT)  is  defined  as: 


+L 

R(kT)=  lim  _i_  x(nT)x(nT+kT ) 

°o  2LT  n  =  -L 


-L  <  k  <  +L 


(2.8) 


where  k  is  the  number  of  samples  by  which  the  time  series 
x(nT )  is  separated  from  its  "lagged"  version  x(nT+kT).  A 
commonly  used  estimator  R(kT)  for  the  autocorrelat ion 
function  of  an  infinite  series  based  on  a  finite  sample  of 
length  N  is: 


N-  Ikl  -1 
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/R(kT )  =  J. 

N 


n=0 


x(nT)x(nT+kT) 


-N+1  <  k  <  N-1 

(2.9) 

It  is  well  known  that  the  power  spectral  density  S(mF) 
and  autocorrelat ion  function  R(kT)  are  related  by  the  DFT 
relat ion: 


/  v  ,  v  -  ]  2TT  no  F  kT 

S(mF)=  7"  R(kT)e  d 


k=- 


OO 


-  °°  <  m  <  +  00 

(2.10) 

Logically  a  suitable  estimator  ^(mF)  for  the  power  spectrum 
of  an  infinite  process  based  on  a  finite  sample  of  length  N, 
is  related  to  the  autocorrelat ion  function  of  the  sample 
R(kT)  by: 


S(mF)  = 


N-1  a  -i2TrmFkT 
ZZ  R(kT)ed 
k=-N+1 


0  <  m  <  N-1 

(2.11) 

If  equation  (2.9)  is  substituted  into  equation  (2.11) 
the  result  may  be  written  in  terms  of  X(mF),  the  discrete 
Fourier  coefficients  of  the  time  series  x(kT)  (ref.  equation 


.  ,  -  -  p,  ■ 
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(2.6)).  The  power  spectrum  estimate  is  then  given  by: 


S(mF) = ( 1/N) | X(mF) | 2 


(2.12) 


The  advent  of  the  FFT  algorithm  has  brought  this  form 
of  the  estimator  into  widespread  usage.  'S(mF),  sometimes 
called  the  periodogram,  indicates  the  distribution  in 
frequency  of  the  power  present  in  the  sample  x(kT).  Hence, 
the  periodogram  will  elicit  the  presence  of  any  strong 
periodic  components.  Care  must  be  taken  in  the 

interpretat ion  of  the  periodogram,  however,  as  there  are  two 

A/  * 

problems  associated  with  the  use  of  S(mF)  as  a  power 
spectrum  estimator: 

(1)  The  time  series  being  analyzed  is  only  a  finite 
sample  of  a  time  series  of  theoretically  infinite  length. 
When  the  power  spectrum  of  the  infinite  time  series  is 
estimated  from  this  sample,  the  finite  length  of  the  sample 
leads  to  the  problem  known  as  leakage. 

(2)  'S(mF)  is  an  inconsistent  estimator  of  the  power 
spectrum  of  a  random  and  infinite  time  series. 

Suitable  techniques  may  be  employed  to  reduce  the 
detrimental  effects  of  each. 

2.4.2  Leakage 

For  practical  reasons,  the  power  spectrum  of  a  time 
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Fiqure  2.3  Effect  of  finite  length  of  data. 
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series  of  theoretically  infinite  length  is  estimated  from  a 
sample  of  finite  length.  In  effect,  the  finite  length  sample 
is  obtained  by  viewing  the  time  series  of  infinite  length 
through  a  rectangular  data  window  of  finite  length  and  unit 
amplitude  (fig.  2.3(e)).  The  continuous  Fourier  transform  of 
the  rectangular  window  has  the  form  of  (sin  x)/x. 
Multiplication  in  the  time  domain  corresponds  to  convolution 
in  the  frequency  domain  (and  vice  versa);  hence  the  Fourier 
transform  of  the  finite  length  of  data  is  equivalent  to  the 
transform  of  the  infinite  length  of  data  convolved  with  the 
transform  of  the  rectangular  window  (fig.  2.3(f)). 
Therefore,  the  calculated  power  at  any  frequency  will  be 
biased  by  leakage  from  adjacent  frequencies.  By  using  a  data 
window  other  than  the  rectangular  ,  it  is  possible  to  better 
1 oca  1 i ze  the  con t r i bu t i on  of  each  f r equency  componen t . 

2.4.3  Inconsistency 

Two  criteria  are  commonly  used  in  the  evaluation  of  an 
estimator:  bias,  and  variance.  Bias  refers  to  the  difference 
between  the  mean  of  a  number  of  separate  estimates  of  the 
value  of  a  certain  statistical  parameter  using  the  same 
estimator,  and  the  expected  value  of  the  parameter.  Variance 
refers  to  the  mean  squared  deviation  of  the  various 
estimates  from  the  expected  parameter  value.  An  unbiased 
estimator  is  one  whose  bias  approaches  zero  as  the  number  of 
estimates  increases;  a  consistent  estimator  is  one  whose 
variance  approaches  zero  as  the  length  of  the  data  used  by 
the  estimator  increases. 


\ 
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The  bias  and  variance  of  spectral  estimates  obtained 
using  equation  (2.12)  have  been  examined.  It  has  been  shown 
(Schwartz  and  Shaw, 1975)  that  the  mean  of  the  distribution 
of  estimates  of  the  autocorrelat ion  estimator  'R(kT)  is 
related  to  R(kT )  by: 


E{R(kT ) } 


( 1-( |k|/N) )R(kT ) 


(2.13) 


where  k  is  the  lag  number  defined  in  equation  (2.8). 

By  making  the  appropriate  substitutions  into  equation 
2.10,  it  can  be  shown  that: 


E{S(mF) }= 


+  oo  -  jSirmFkT 

51  R(kT)w(kT)e 0 
k =-<*> 


where  w(kT)  is  defined  as: 


(2.14) 


k  <  N 
k  >  N 

(2.15) 

It  can  be  seen  by  comparing  equations  (2.9)  and  (2.14) 
that  E{S(mF)}  differs  from  S(mF)  by  the  factor  w(kT). 
Therefore,  the  periodogram  is  a  biased  estimator  of  the 
power  spectrum.  The  term  w(kT)  is  known  as  a  window  function 
and  results  from  the  use  of  data  of  less  than  infinite 
duration.  However,  as  N  is  increased  the  bias  is  reduced. 


w(kT) =  1-( |k|/N) 
w(kT ) =  0 


' 
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Unfortunately,  the  same  cannot  be  said  for  the  variance  of 
the  estimator. 

If  the  data  comes  from  a  Gaussian  process,  it  can  be 
shown  (Jenkins  and  Watts, 1968)  that 

A . 

Lim  var  S(mF)^S2(mF) 

oo 

(2.16) 

i.e.  the  variance  of  S(mF)  approaches  the  square  of  the 
true  spectrum  at  each  frequency  mF.  Equation  (2.16),  which 
applies  approximately  to  non-Gaussian  data  as  well,  implies 
that  calculations  using  different  sets  of  N  samples  from  the 
same  signal  could  yield  radically  different  estimates  of  the 
true  power  spectrum,  even  for  large  N. 

It  can  be  shown  that  the  variance  of  an  estimate  of  the 

\ 

power  spectrum  using  equation  (2.12)  using  N  data  samples 
can  be  reduced  by  segmenting  the  data  into  shorter  pieces 
and  then  averaging  the  periodograms  calculated  from  the 
pieces  (Welch, 1967) .  An  alternate  method  involves  a 
smoothing  operation,  performed  on  the  periodogram  obtained 
from  a  single  complete  record  (Schwartz  and  Shaw, 1975). 

2.5  Convolut ion 

Since  the  multiplication  of  the  Fourier  transforms  of 
two  signals  is  equivalent  to  the  Fourier  transform  of  the 
convolution  of  the  two  signals,  the  process  of  convolution 
can  be  achieved  by  three  Fourier  transformat  ions  and  a 
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Figure  2.4(a)  Aperiodic  series  x(kT)  and  h(kT),  and  their 
aperiodic  convolution  c(kT). 


Figure  2.4(b)  Periodic  series  x(kT)  and  h(kT),  and  their 
periodic  convolution  c(kT).  Intervals  of  interference  are  shaded. 
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multiplication.  For  sufficiently  large  N  it  has  been  shown 
( Stockham, 1966 )  that  convolution  achieved  through  frequency 
transf ormat ion  with  the  FFT  algorithm  requires  much  less 
computing  time  than  performing  the  convolution  in  the  time 
domain.  Therefore,  evaluating  the  convolution  integral  by 
utilizing  the  FFT  algorithm  is  known  as  the  "fast 
convolution"  method. 

A  problem  arises,  however,  since  the  FFT  applies  to  the 
DFT  only,  and  multiplication  of  DFTs  corresponds  to 
convolution  of  periodic  functions,  not  aperiodic  ones.  This 
problem  may  be  overcome  by  padding  both  functions  with  a 
sufficient  number  of  zero  values  to  render  the  periodic 
convolution  into  an  aperiodic  one.  If  x(kT)  is  defined  for 
0  <  k  <  M-1  and  h(kT)  is  defined  for  0  <  k  <  N-1,  then  each 
of  the  functions  must  be  padded  out  to  a  new  length  L,  where 
L  >  M+N-1,  to  render  the  periodic  convolution  into  an 
aperiodic  one  (fig.  2.4).  It  is  also  necessary  for  L  to  be 
an  integer  power  of  two  to  obtain  the  maximum  computational 
expediency  possible  with  the  FFT  algorithm. 


2.6  F i Iter  design 

2.6.1  Introduction 

In  order  to  obtain  estimates  of  the  particles7 
velocities,  it  was  necessary  to  differentiate  the 
position-time  data.  As  a  conventional  differentiation 
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process  was  expected  to  cause  undue  emphasis  of  the  random 
high  frequency  components  present  in  the  time  series 
representing  particle  trajector ies ,  the  need  for  low  pass 
filtering  was  anticipated.  The  design  of  abandlimited 
di f ferent iator  was  therefore  investigated. 

As  shown  in  equation  (2.2),  the  inverse  Fourier 
transform  of  X(f)  is  given  by: 

X(f)  ed  df 

The  derivative  of  x(t)  is  therefore  given  by: 


•+ 

x(t)=  X 


d 

dt 


4-  CO 


I  2TT -Tt 

{  x ( t )  }  =  d  f  X(f)  e 4  df  } 
dt  J 


-OO 


4-c»o 


•  r  i2ir-ft 

=  X(f)  die6  }  df 
J  dt 


-  oo 

+  CO 


T  iCTft 

=  J  ( j2TT f )  X(f)  e1*  df 


—  oO 


(2.17) 


Thus  it  can  be  seen  that  differentiation  with  respect  to 
time  of  a  function  of  time  x(t)  may  be  performed  in  the 
frequency  domain  by  multiplying  its  Fourier  transform  X(f) 
by  the  function  H(f),  where  H(f)=  j2TTf,  and  then  computing 
the  inverse  Fourier  transform  of  the  product.  The  expression 
for  a  bandlimited  differentiator  in  the  frequency  domain  is 
given  by: 
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H(  f  )  =  j2TTf 
=  0 


| f |  <  Fcut 
otherwise 

(2.18) 


where  Fcut  is  the  filter  cutoff  frequency. 

For  a  number  of  practical  reasons  ideal  filters  like 
the  above  can  only  be  approximated.  Methods  such  as 
frequency  sampling,  ripple  equalization,  and  Fourier 
expansion  and  windowing  have  all  been  used  (Childers  and 
Dur 1 ing, 1975;  Antoniou, 1979) .  Each  method  uses  a  different 
criterion  to  define  the  "best"  approximation  to  a  given 
"ideal"  filter.  They  all,  however,  lend  themselves  to 
application  in  either  the  time  or  the  frequency  domain,  with 
identical  results. 

2.6.2  Fi Iter ing  with  fast  convolut ion 

Implicit  in  the  principles  of  fast  convolution  is  the 
periodicity  of  the  functions  being  convolved,  and  hence  the 
result  of  fast  convolution  is  periodic,  rather  than 
aperiodic  convolution.  Unless  both  functions  are  of  finite 
duration  the  condition  L  >  M+N-1  discussed  in  Section  2.5  as 
being  required  to  render  the  periodic  convolution  into  an 
aperiodic  equivalent,  cannot  be  met  for  any  L. 

On  initial  examination  it  would  appear  that  "ideal" 
filtering  could  be  performed  in  the  frequency  domain.  One 
need  only  specify  the  continuous  frequency  response  of  the 
desired  filter,  sample  this  response  to  obtain  the  Fourier 
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(a)  Frequency  domain  coefficients  of  la-.'  pass  differentiator 


(b)  Impulse  response  of  filter  specified  in  (a) 


Figure  2 . 5 
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coefficients  of  the  discrete  version  of  the  filter  in  the 
frequency  domain,  multiply  the  Fourier  transform  of  the 
signal  to  be  filtered  by  these  coefficients,  and  inverse 
transform  the  product  to  obtain  the  filtered  output  in  the 
time  domain. 

If,  however,  one  computes  the  inverse  Fourier  transform 
of  the  sampled  frequency  response  of  an  ideal  low  pass 
differentiating  filter,  a  periodic  impulse  response  is 
obtained.  The  periodic  nature  of  this  response  results  in  an 
effect  which  is  similar  to  the  aliasing  problem  in  the 
frequency  domain  caused  by  undersampling  a  signal  in  the 
time  domain.  The  tails  of  each  impulse  response  interfere 
with  those  of  adjoining  impulse  responses  due  to  their 
infinite  length.  Figure  2.5(a)  illustrates  the  imaginary 
components  of  the  frequency  response  of  a  low  pass 
dif ferent iator  specified  in  equation  (2.18)  with  a  cutoff 
frequency  of  .3  Hz,  whose  frequency  response  has  been 
sampled  NT  times  in  the  range  from  from  0  to  1/T  Hz  at 
intervals  of  1/NtT  Hz.  (Note  that  for  the  case  of  the 
differentiating  filter  all  the  real  coefficients  in  the 
frequency  domain  are  zero).  Figure  2.5(b)  illustrates  the 
periodic  sampled  impulse  response  of  this  low  pass 
different iator ,  and  was  obtained  by  computing  the  inverse 
DFT  of  the  NT  samples  of  the  frequency  response.  The 
interference  of  the  adjoining  responses  is  apparent. 

This  interference  can  only  be  eliminated  by  truncating 
the  impulse  response  of  the  filter  (in  effect,  cutting  off 
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the  tails).  Truncation,  however,  gives  rise  to  Gibb's 
phenomenon,  which  manifests  itself  as  a  fixed  percentage 
overshoot  and  ripple  before  and  after  any  sudden  transitions 
in  the  frequency  response.  In  the  case  of  the  low  pass 
dif ferent iator ,  this  occurs  at  the  cutoff  frequency,  where 
the  function  H(f)  abruptly  goes  from  H(f)=j27Tf  to  H(f)=0. 

Considerable  efforts  have  been  made  by  workers  in  the 
field  to  decrease  the  overshoot  and  ripple.  It  has  been 
shown  that  "windowing"  the  truncated  inpulse  response  can 
sign i f icant ly  reduce  the  effect  of  Gibb's  phenomenon  (Kuo 
and  Kaiser, 1966;  Helms, 1968;  Rabiner  et  al.,1970).  The 
truncated  and  windowed  inpulse  response  can  sometimes  be 
made  to  yield  an  acceptable  approximation  to  the  "ideal" 
filter  originally  specified. 

The  shape  of  window  used  will  obviously  alter  the 
filter  character ist ics.  Selection  of  an  appropriate  window 
is  usually  based  on  consideration  of  the  trade-offs  inherent 
in  minimizing  the  duration  of  the  filter's  impulse  response, 
the  ripple  in  the  filter's  frequency  response,  and  the  width 
of  the  region  of  filter  response  lying  between  passband  and 
reject  band  known  as  the  transition  width.  As  windows  are  a 
major  topic  in  Fourier  analysis,  their  use  was  examined  in 
some  detai 1 . 

2.6.3  W i ndow  carpentry 

The  topic  of  window  design  has  been  widely  explored  in 
the  literature  (Childers  and  Dur 1 ing, 1975;  Bloomf ield, 1976 ; 
Hamming, 1977) .  The  objective  is  usually  to  determine  a  time 
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(a)  Tukey  and  rectangular  data  windows 
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(b)  Carparison  of  the  frequency  response  of  the  Tukey  and  rectangular  windows 


Figure  2.6  The  Tukey  and  rectangular  data  windows 
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domain  window  which  has  lower  side lobes  in  the  frequency 
domain  than  does  the  rectangular  window.  Windows  could  be 
used  at  several  points  in  the  analytical  procedure: 

(1)  The  "Tukey"  data  window,  among  others, is  often  used 
in  power  spectrum  estimation  to  reduce  the  effect  of  leakage 
through  the  sidelobes  of  the  rectangular  data  window  in  the 
frequency  domain.  This  window  consists  of  a  raised  cosine 
shaped  taper  of  the  first  and  last  10%  of  the  data,  and  a 
constant  value  of  one  in  between.  The  smaller  sidelobes  of 
this  data  window  in  the  frequency  domain  when  compared  with 
those  of  the  rectangular  window  are  obtained  at  the  cost  of 
a  slight  increase  in  the  main  lobe  bandwidth.  The  continuous 
time  functions  of  these  two  windows  were  each  sampled  NF 
times  at  T  second  intervals  to  produce  the  discrete  time 
funct ions  shown  in  figure  2.6(a).  It  should  be  noted  that 
the  frequency  responses  of  discrete  time  windows  are 
continuous  and  periodic.  The  frequency  responses  of  these 
two  windows  are  illustrated  in  figure  2.6(b). 

(2)  A  window  can  be  used  to  modify  the  impulse  response 
of  a  filter  in  the  filter  design  process.  The  Kaiser  window 
has  been  shown  to  have  character ist ics  well  suited  to  this 
application.  By  adjusting  oc  ,  the  shape  parameter  of  the 
window,  the  sidelobes  of  the  window's  frequency  response  can 
be  diminished  at  the  cost  of  increased  main  lobe  bandwidth. 

2.6.4  Ka i ser  w i ndow 

An  ideal  window  would  be  one  which,  in  the  frequency 
domain,  had  a  main  lobe  of  negligible  width,  and  sidelobes 
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(a)  The  effect  of  the  shape  parameter  on  the  Kaiser  window 
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(b)  The  effect  of  the  shape  parameter  on  the  frequency  response  of  the  Kaiser  window 


Figure  2.7  The  Kaiser  windcw 
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of  negligible  magnitude.  In  practice,  this  is  impossible  to 
achieve,  as  reduction  of  main  lobe  width  can  be  obtained 
only  at  a  cost  of  increased  side lobe  magnitude.  However, 
J.F.  Kaiser  has  developed  a  flexible  family  of  window 
functions  whose  frequency  responses  have  a  narrower  mainlobe 
bandwidth  for  a  given  magnitude  of  side lobes  than  any  other 
window  function  presently  known  (Kuo  and  Kaiser, 1966; 
Hamming, 1977 ) .  This  function,  called  the  Kaiser  window,  is 
given  in  discrete  form  by: 


w(kT)=  !»{=<■  Vi  -  ( 2k/  (Nf-  1 ) ) 2  } 

=  0 


| k |  S(Np-1)/2 


otherwise 


(2.19) 


where  I0  is  the  modified  Bessel  function  of  the  first  kind 
of  order  zero,  Np  is  the  number  of  terms  comprising  the 
window,  and  the  parameter  o c  is  the  shape  factor.  The  width 
of  the  Kaiser  window  is  denoted  by  D,  and  is  related  to  Np 
and  T,  the  period  between  the  terms,  by  the  equation: 


D=  ( Np - 1 ) T 


(2.20) 


The  frequency  response  of  w(kT),  W(f),  is  given  by: 
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W(f)=  D  sinh  {oC  VT- 


2 


=  D  sin  {o C  /(f/fa  )z-l  } 

P-1 


D 

ITTScT 


f  <  fq 


f  >  fq 

f  =  fo 


(2.21) 


where 


f*  =  <*/TTD. 


(2.22) 


The  effect  of  varying  the  shape  parameter  on  a  Kaiser  window 
of  width  D  seconds  sampled  Np  times  at  intervals  of  T 
seconds,  and  on  its  frequency  response  may  be  seen  in  figure 
2.7. 

2.7  Digital  f i 1  ter  design 

Several  parameters  are  usually  important  in  the  design 
of  digital  filters  using  the  Kaiser  window:  D,  the  duration 
of  the  filter's  impulse  response;  Af,  the  transition  width 
between  the  reject  and  pass  band  in  the  filter's  frequency 
response;  and  5  ,  the  magnitude  of  the  maximum  ripple  in  the 
filter's  frequency  response  normalized  to  the  maximum 
magnitude  of  the  response.  The  three  parameters  are  not 
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independent,  making  it  necessary  to  specify  only  two  of  the 
three.  The  parameters  are  related  by  the  following  equation: 


D=  2( A-7.95) T 
28.72  Af 


(2.23) 


where  A  is  known  as  the  attenuation  factor,  and  gives  the 
magnitude  of  the  maximum  response  in  the  filter's  frequency 
response  relative  to  that  of  the  maximum  ripple,  in  units  of 
dB.  A  is  related  to  S  by: 


A=  -20  log  (  S  ) 


(2.24) 


A  formula  has  been  empirically  determined  by  Kaiser 
(Kaiser , 1974;Hamming, 1977)  relating  the  window  shape 
parameter  o t  to  the  attenuation  factor  A. 


=  . 1 102( A-8. 7 )  A>50 

=  .  5842(  A-21  )*4  +  .  07886(  A-21 )  21<A<50 

=  0  A<21 

(2.25) 


These  equations  may  be  used  in  conjunction  with  the 
Fourier  Transform  and  the  Kaiser  window  to  perform  filter 
design.  The  method  may  be  easily  implemented  on  a  digital 
computer,  and  despite  its  computational  simplicity,  the 
resultant  nonrecursive  filter  compares  very  favourably  with 
filters  designed  with  much  more  involved  calculations,  such 
as  the  least  squares,  Chebyshev,  and  other  methods  (Kuo  and 
Kaiser, 1966;  Rabiner  et  a  1., 1970;  Rabiner , 1971 ) . 
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Assume  it  is  necessary  to  filter  a  time  series 
consisting  of  M  equally  spaced  data  points  sampled  at  T 
second  intervals,  utilizing  fast  convolution.  Using  the 
Kaiser  window  and  the  Discrete  Fourier  Transform  as  tools,  a 
suitable  filter  would  be  designed  as  follows  (refer  to 
f igure  2.8) : 

(a)  two  of  the  three  filter  character ist ics-  transition 
width  Af,  maximum  normalized  ripple  &  in  the  frequency 
response,  and  duration  D  of  the  impulse  response-  are 
specified.  The  third  parameter  is  calculated  from  the  other 
two,  using  equations  2.22  and  2.23.  The  window  shape 
parameter  ©6  is  calculated  from  the  maximum  ripple 
specification,  using  Eq.  2.24.  D  is  constrained  by  the 
condition  necessary  to  render  the  periodic  convolution  into 
an  aperiodic  equivalent  that  D  <  (L-M)T,  where  L  is  the 
number  of  points  in  the  Fourier  transform  to  be  used  in 
performing  the  filtering  via  fast  convolution. 

(b) ,  the  frequency  response  of  the  desired  filter  is 
specified  over  the  range  0  Hz  to  1/T  Hz,  and  is  sampled  Ny 
times  at  intervals  of  1/NyT  Hz,  where  Ny  is  the  number  of 
points  in  the  Fourier  transform  to  be  used  in  the  filter 
design  process.  If  the  impulse  response  of  the  resultant 
filter  is  to  be  limited  to  D  seconds,  it  is  suggested  that 
Ny  >  1 6 ( D/T  +1)  (Rabiner  et  al.,1970).  This  will  reduce  the 
error  in  estimating  the  filter  impulse  response  h(kT)  caused 
by  the  interference  of  adjacent  tails  of  the  periodic 
impulse  response. 
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(c)  ,  the  inverse  DFT  of  the  NT  samples  of  the  frequency 
response  is  computed,  resulting  in  the  time  domain  samples 
h(kT)  of  the  filter  impulse  response.  The  impulse  response 
is  then  truncated  to  a  duration  of  D. 

(d)  ,  the  samples  w(kT)  of  the  Kaiser  window  of  width  D 
seconds  are  calculated  using  equation  2.18,  with  a  period 
between  terms  of  T  seconds.  The  truncated  filter  impulse 
response  samples  h(kT)  computed  in  (c)  are  tapered  with  this 
w i ndow . 

(e) ,  the  result  is  h#(kT),  the  sampled  impulse  response 
of  the  specified  filter.  The  impulse  response  h'  (kT)  is 
padded  with  zero  terms  to  a  total  length  of  L  points,  and 
the  filtered  output  is  then  obtained  by  fast  convolving  the 
data  of  L  points  with  the  filter  impulse  response  of  L 
points. 
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CHAPTER  III 


FILTER  DESIGN  AND  EVALUATION 

3. 1  Introduction 

A  band  limited  differentiating  filter  was  designed 
using  the  windowed  Fourier  expansion  method.  Simulated 
particle  trajectories  were  analyzed  with  this  filter,  and 
the  results  thus  obtained  compared  with  those  obtained  using 
a  two  term  approximation  to  an  ideal  di f ferent iator . 

3.2  Design 

D,  the  duration  of  the  impulse  response  of  the  filter, 
was  limited  by  the  formula  D  <  (L-M)T,  the  condition 
necessary  to  render  the  periodic  convolution  into  an 
aperiodic  equivalent.  The  longest  trajectory  data  set 
consisted  of  350  samples  (T=  1/3  second);  therefore  M,  the 


number  of  samples  in 

the  longest 

data 

set 

to 

be  filtered, 

was  conser va  t i ve 1 y 

chosen  to 

be 

400. 

L, 

the  number  of 

Fourier  coefficients 

to  be  used 

in 

the 

fast 

convolut ion 

procedure,  was  limited  by  available  memory  space  to  1024 
terms,  allowing  a  frequency  resolution  of  about  .003  Hz.  As 
a  result,  D  was  limited  to  a  maximum  of  208  seconds 
(Np=625).  Therefore,  an  impulse  response  of  up  to  625 
samples  could  be  convolved  with  a  data  set  of  up  to  400 
samples,  with  no  error  introduced  by  cyclic  convolution. 

In  the  majority  of  applications,  a  filter  with  a 
relatively  short  impulse  response  is  used  to  process  a 
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signal  of  comparatively  long  duration.  In  this  way  the 
errors  due  to  the  so-called  "end  effects"  appear  in  only  a 
small  percentage  of  the  output.  However,  a  reduction  in  the 
number  of  terms  comprising  the  filter's  impulse  response 
results  in  a  less  accurate  frequency  response  when  compared 
with  that  of  an  ideal  filter.  Errors  due  to  end  effects  were 
reduced  by  shortening  the  filter  impulse  response,  at  the 
expense  of  increasing  the  error  of  the  filter  frequency 
response.  A  filter  with  a  impulse  response  of  Np  terms  will 
introduce  end  effect  errors  into  Np-1  terms  of  the  output. 
Since  the  data  being  filtered  in  this  application  was  of 
short  length  (34-350  readings),  the  filter  used  had  an 
impulse  response  of  much  less  than  625  terms.  Since 
eighty-eight  of  the  one  hundred  seventy-two  trajectory  data 
sets  contained  more  than  100  data  points,  it  was  decided  to 
design  the  filter  such  that  end  effects  appeared  in  a 
maximum  of  20  terms  of  the  output.  (Thus,  end  effects 
appeared  in  no  more  than  20%  of  the  results  obtained  from 
each  of  eighty-eight  trajectories).  A  filter  with  an  impulse 
response  of  <  21  terms  was  required. 

The  filter  design  procedures  utilizing  the  Kaiser 
windowed  Fourier  expansion  method  were  implemented  in  HP 
FORTRAN  IV  on  a  Hewlett  Packard  Model  2100A  minicomputer. 
Available  memory  space  limited  the  length  of  the  Fourier 
transform  used  in  the  filter  design  procedures  to  a  maximum 
of  2048  terms  thus  determining  Nj.  Rabiner's  suggestion  that 
NTbe  chosen  >  16Np  to  keep  aliasing  error  in  the  design 
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procedure  to  an  insignificant  level  (Rabiner  et  al.,  1970), 
was  therefore  easily  met  in  all  cases.  Filters  composed  of 
21,  31,  and  41  terms  were  designed.  The  remaining  window 
parameters  were  calculated  from  equations  (2.23),  (2.24), 
and  (2.25).  Some  possible  combinations  are  shown  below: 


A(db) 

Np 

Af (Hz) 

20 

21 

.042 

0.0 

20 

31 

.028 

0.0 

20 

41 

.021 

0.0 

30 

21 

.077 

2. 12 

30 

31 

.051 

2.12 

30 

41 

.038 

2.12 

40 

21 

.112 

3.4 

40 

31 

.074 

3.4 

40 

41 

.056 

3.4 

Table  3.1 


The  nine  filters  specified  in  the  preceding  table  were 
designed  to  perform  band  limited  differentiation  on  data 
sampled  at  3  Hz.  The  cutoff  frequency  was  chosen  to  be  .3 
Hz,  based  on  consideration  of  the  signals7  power  spectra 
(ref.  Section  5.3).  The  frequency  responses  of  the  filters 
specified  in  table  3.1  are  plotted  in  figure  3.1.  Note  that 
although  the  frequency  responses  are  periodic  only  the 
region  of  interest  (i.e.  from  0  to  1.5  Hz)  has  been  plotted. 
The  trade-offs  among  reduction  of  the  transition  width  and 
ripple  in  the  frequency  response,  and  of  the  duration  of  the 
impulse  reponse  in  the  time  domain,  are  apparent. 

The  "best"  filter  to  use  in  the  analysis  of  particle 
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(a)  A=  20  dB. 


I  €  .  0 


(b) 


(Hz) 

A=  30  dB. 


(c)  A=  40  dB. 


Figure  3.1  Corrparison  of  filter  frequency  responses 
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motion  was  not  cleanly  indicated.  Since  the  cutoff  frequency 
used  in  the  low  pass  filtering  was  somewhat  arbitrarily 
chosen,  the  transition  width  was  not  critical.  This  would 
not  have  been  the  case  had  it  been  necessary  to  eliminate 
strong  components  falling  only  slightly  above  the  cutoff 
frequency.  In  all  of  the  cases  examined,  at  the  cutoff 
frequency  the  trajectory  power  spectrum  had  dropped  several 
orders  of  magnitude  from  its  peak  power,  hence  the 
transition  width  was  of  minor  importance.  Three  of  the  nine 
designs  were  selected  to  be  evaluated  through  the  use  of 
simulations.  Filters  whose  maximum  normalized  ripples  were 
approximately  30  dB  below  the  magnitude  of  the  maximum 
response  in  the  frequency  response  were  chosen. 

3 . 3  Simulat ions 

Simulated  trajectories  generated  by  a  summation  of  sine 
and  cosine  terms  to  which  an  average  velocity  component  had 
been  added  closely  resembled  actual  particle  trajectories, 
and  were  used  to  estimate  the  error  inherent  in  the 
instantaneous  velocity  estimation  methods.  Simulation 
parameters  were  based  on  observations  of  real  particles7  x 
and  y  average  velocities  (Vx  and  V y  respectively),  and  on 
the  frequency  and  magnitude  of  periodic  components  (as 
evidenced  by  the  power  spectrum) . 

For  example,  the  position  coordinates  at  any  time  t 


could  be  given  by: 


' 


eft 
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x(t)=  2.0  cos(.05t)  +  .5  sin(.02t)  +  Vx  t 
y(t)=  .1  cos(.03t)  +  .08  sin(.025t)  +  t 
d(t)=  /{x( t )-x(0)}2+{y(t)-y(0)}2 

and  the  velocities,  through  differentiation,  were  given  by 

Vx(t)=  -.1  sin(.05t)  +  .01  cos(.02t)  +  Vx 
Vy ( t ) =  -.03  sin(.03t)  +  .002  cos(.025t)  +  Vy 

Vd ( t )  =  vyftwm  +  \i\,i  tiwrti 

/{x( t ) -x( 0) } 2  +  {y ( t ) -y ( 0 ) } 2 

A  finite  number  of  samples  N  were  taken  at  T  second 
intervals  from  each  of  the  functions  x(t),  y(t),  and  d(t). 
In  some  cases  the  effect  of  noise  was  simulated  by  adding  to 
each  sample  a  random  number  from  a  Gaussian  distribution  of 
numbers  with  a  specified  mean  and  standard  deviation.  With  a 
distribution  mean  of  zero,  this  population  contains  equal 
power  at  all  frequencies.  The  standard  deviation  of  the 
population  was  made  sufficiently  large  to  ensure  that  the 
power  present  above  .3  Hz  in  the  new  “noisy"  trajectory  was 
at  least  as  great  as  that  seen  in  a  real  trajectory.  The  new 
noisy  displacement  samples  were  computed  from  the  noisy  x 
and  y  position  samples.  The  instantaneous  velocity 
components  were  estimated  from  these  samples  using  a  variety 
of  methods.  The  estimates  thus  obtained  were  compared  with 
the  true  velocities  as  given  by  the  expressions  derived  for 
Vx,  Vy,  and  Vd.  The  error  in  the  estimates  was  calculated 
from  the  difference  between  the  estimated  and  the  true 


. 
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velocity  at  each  of  the  N  sample  points.  The  mean  and 
standard  deviation  of  the  N  differences  were  calculated  in 
each  case.  As  the  standard  deviation  of  the  differences  was 
consistently  at  least  one  order  of  magnitude  larger  than  the 
mean,  it  was  chosen  as  the  measure  of  error  in  the  method. 
The  error  of  a  method  shall  henceforth  refer  to  the  standard 
deviation  of  the  error  in  the  estimates  it  produces.  The 
objective  was  to  select  the  method  with  the  smallest  error. 

Figure  3.2  illustrates  the  posi t ion- t ime  data  for  an 
actual  particle.  Figure  3.3  illustrates  this  data  for  a 
simulated  trajectory,  and  figure  3.4  illustrates  this  data 
for  the  simulated  trajectory  with  added  noise.  In  figure  3.5 
the  estimated  power  spectra  of  the  sampled  sets  x(kT)  from 
the  three  trajectories  after  trend  removal  and  Tukey 
windowing  may  be  compared.  Note  that  for  increased 
resolution  the  portions  of  the  spectra  above  .3  Hz  have  been 
plotted  on  a  separate  scale. 

3.4  F i 1  ter  and  procedure  select  ion 

The  velocities  of  a  number  of  simulated  particles  were 
estimated  with  a  variety  of  methods  in  order  to  select  the 
best  of  these  for  application  to  real  particles.  The  results 
for  a  typical  series  of  analyses  are  presented  in  tables  3.2 
and  3.3.  Using  minimization  of  the  method  error  as  the 
criterion,  filter  'A'  and  procedure  (b)  were  selected  as  the 
best  method. 


' 
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I.  F  i  1  tens 
FILTER  A 

Np  =21 

A=30  dB. 

A  f = . 077  Hz 
©C  =2.12 


FILTER  B 
Np=31 
A=30  dB. 

A  f = . 051  Hz 
oC  =2.12 


FILTER  C 
Np=41 
A=30  dB. 

A  f = . 038  Hz 
oC  =2.12 


1 1  F i 1  ter ina  procedures 

(a)  remove  trend  from  data,  convolve  modified  data  with 
f i Iter . 

(b)  remove  trend  from  data,  convolve  modified  data  with 
filter,  discard  first  and  last  10  terms  of  the  results. 

(c)  remove  trend  from  data,  cosine  taper  first  and  last  10 
terms  of  the  data,  convolve  modified  data  with  filter. 

(d)  remove  trend  from  data,  cosine  taper  first  and  last  10 
terms  of  the  data,  convolve  modified  data  with  filter, 
discard  first  and  last  10  terms  of  the  results. 


Ill  Simulat ions 

Case  3.4(a) 

x ( t )  =  Vx  t  +  Bx  +  ll  {Amcos(27T  f  l^t )  +  B^sin(2TT  f  1mt ) } 

I 

y ( t )  =  Vy  t  +  By  +  21  {CnCos( 21T f 20t )  +  Dnsin(2TT f2nt ) } 

n-  » _ 

d ( t ) =  /{x( t ) -x(0) } 2  +  {y(t)-y(0)}2 

Vx = 1 . 0  Vy = . 04 
Bx=0.0  By=5.0 


A, =1.0 
A2=.8 

03  CD 

ro  — 

ii  ii 

f 1,  =.03 
f 12=-05 

C»  =  .2 

D 1  =  .  3 

f 2  j  =  . 02 

C2=.3 

Dz  =  .5 

f 2o= • 08 

A3=.4 
A4  =  .6 
A5  =  •  2 

B3=.3 
B4  = .  5 
B5=.  1 

f 13=.08 

f if  =  .  1 
f  I5  =  ■  13 

C3=.4 

D3=.3 

f23=.125 

Positions  sampled  N=100  times  at  interval  T=1/3  second 

case  3.4(b) 

x7  (kT)=xIkT)+N(0-0,0.13) 
y7  (kT)=y(kT)+N( 0.0,0. 13) 

d'  (kT)=y(x'  (kT)  -x#  (0) ) 2+(y'  (kT) -y7  (0) ) 2 

N(a,b)-  function  that  generates  pseudo-random  Gaussian  noise 
of  mean  ' a7  and  standard  deviation  7 b7  . 


.  . 
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Case  3.4(a)  results 


Filter 

Est imate 

method 

Vx  error 
(um/sec) 

V y  error 
(um/sec) 

Vd  error 
(um/sec) 

Mean 

S.D. 

Mean 

S.D. 

Mean 

S.D. 

A 

(a) 

-2.9E-2 

1 .3E-1 

+6.3E-4 

2.8E-2 

-3.0E-2 

1 . 4E- 1 

A 

(b) 

+7.7E-4 

1 .3E-2 

+9.5E-4 

1.3E-2 

+1.1E-3 

1 . 4E-2 

A 

(c) 

-6. IE-2 

2.2E-1 

-4.7E-3 

6.2E-2 

-6. IE-2 

2.2E-1 

A 

(d) 

+2.7E-3 

2.4E-2 

+8.8E-4 

1 .3E-2 

+3.0E-3 

2.5E-2 

B 

(a) 

-3.0E-2 

1 . 4E- 1 

-6.9E-5 

2.7E-2 

-3.0E-2 

1 . 4E- 1 

B 

(b) 

-3.6E-3 

1 .6E-2 

-3.7E-4 

6.4E-3 

-3.4E-3 

1 .5E-2 

B 

(c) 

-6.2E-2 

2.3E-1 

-5.4E-3 

6.4E-2 

-6.2E-2 

2.2E-1 

B 

(d) 

+2.4E-4 

2. IE-2 

-8.2E-6 

8.8E-3 

+4.4E-4 

2.3E-2 

C 

(a) 

-3.0E-2 

1 .4E-1 

+5.6E-4 

2.7E-2 

-3.0E-2 

1 .5E-1 

C 

(b) 

-2.4E-3 

2.4E-2 

-1 .4E-5 

7.2E-3 

-2.2E-3 

2.2E-2 

C 

(c) 

-6.2E-2 

2.3E-1 

-4.7E-3 

6.4E-2 

-6.2E-2 

2.3E-1 

C 

(d) 

+2.5E-3 

2.0E-2 

+5.7E-4 

8.7E-3 

+2.6E-3 

2.2E-2 

Table 

3.2 

Case 

3.4(b) 

results 

F i 1  ter 

Est imate 

Vx  error 
(um/sec) 

Vy  error 
(um/sec) 

Vd  error 
(um/sec) 

method 

Mean 

S.D. 

Mean 

S.D. 

Mean 

S.D. 

A 

(b) 

-1 .4E-3 

5.0E-2 

-1.1E-3 

5.6E-2 

-1 . IE-3 

4.7E-2 

A 

(d) 

+8.6E-4 

5.9E-2 

-6.5E-4 

5.7E-2 

+1.1E-3 

5.6E-2 

B 

(b) 

-6.0E-3 

5.3E-2 

-2.2E-3 

5.3E-2 

-5.8E-3 

5.2E-2 

B 

(d) 

-1 .6E-3 

6.0E-2 

-1 .3E-3 

5.4E-2 

-1 .5E-3 

5.8E-2 

C 

(b) 

-4.7E-3 

5.7E-2 

-1 .9E-3 

5.6E-2 

-4.6E-3 

5.6E-2 

C 

(d) 

+7.6E-4 

5.9E-2 

-6.7E-4 

5.6E-2 

+8.4E-4 

5.8E-2 

Table  3.3 
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0.0  8.0  16.0  24.0  32.0 

POSITION  LUtt) 

(a)  Trajectory  S4 


(b)  X  position  versus  time 


(d)  X  deviations  from  straight  line  fit 


O 
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(c)  Y  position  versus  time 


CD 

LD 


(e)  Y  deviations  from  straight  line  fit 


Figure  3.2  Position-time  data  for  a  real  trajectory. 
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(b)  X  position  versus  time 
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(c)  Y  position  versus  time 
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TIME  t  SEC  J 

(d)  X  deviations  from  straight  line  fit 
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Figure  3.3  Position-time  data  for  a  simulated  trajectory 
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(b)  X  position  versus  time 


□ 

CD 


(c)  Y  position  versus  time 


□ 

CD 


.  U  L.  U  • 

TIME  l SEC ) 


(e)  Y  deviations  from  straight  line  fit 


Figure  3.4  Position-time  data  for  simulated  trajectory  with  noise 
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(a>  Power  spectruir.  cf  actual  trajector- 


Figure  3.5  Power  spectra  of  the  x  component 
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3 . 5  Reduct  ion  of  end  effect  errors 

The  error  in  a  velocity  estimation  method  due  to  filter 
end  effects  can  be  reduced  by  employing  a  filter  whose 
impulse  response  is  of  a  comparat ively  short  duration 
compared  to  that  of  the  data  being  filtered.  It  is  also 
possible  to  reduce  the  error  due  to  end  effects  by 
discarding  a  number  of  the  end  points  from  the  velocity 
estimates  of  a  particle. 

The  following  procedure  was  used  to  assess  end  effects. 
Estimates  of  Vx,  Vy,  and  Vd  were  obtained  from  100  samples 
(the  typical  number  of  samples  obtained  from  an  actual 
trajectory)  using  the  method  previously  described.  After 
discarding  a  given  number  of  result  endpoints,  the  method 
error  was  calculated  from  the  remaining  estimates  as  before. 
Figure  3.6  illustrates  the  effect  of  varying  the  number  of 
endpoints  discarded  on  the  method  error  in  Vx,  with  and 
without  noise  added  to  the  trajectory  samples.  Clearly, 
discarding  up  to  five  points  from  each  end  of  the  estimates 
produced  a  substantial  reduction  of  the  error.  However, 
discarding  more  than  five  points  from  each  end  achieved 
little  further  reduction. 

3.6  Velocity  est imat ion  procedure 

Based  on  the  results  from  the  preceding  trials,  the 
following  method  (later  referred  to  as  'method  2' )  was  used 
to  estimate  particle  instantaneous  velocity  components: 

(a)  the  trend  as  estimated  by  the  linear  regression 
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Figure  3.6  Effect  of  discarding  result  endpoints  on 
the  error  in  the  estimates  of  Vx . 
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line  was  subtracted  from  the  samples  of  the  desired  position 
component  data. 

(b)  the  trend-free  samples  were  padded  with  zeros  to 
1024  points,  and  fast  convolved  with  a  21  term  band  limited 
differentiating  filter  (filter  7 A#  ) ,  using  a  Fourier 
transform  of  1024  points. 

(c)  five  points  were  discarded  from  each  end  of  the 
velocity  estimates  obtained  in  (b). 

(d)  the  average  velocity  estimate  given  by  the  slope  of 
the  linear  regression  line  in  (a)  was  added  to  the  estimates 
obtained  in  (c).  The  results  were  the  estimates  of  the 
selected  component  of  the  particle's  instantaneous  velocity. 


3.7  Comparison  with  the  two  term  approximation 

After  establishing  the  method  by  which  velocity 
estimation  would  be  undertaken,  it  was  desirable  to  compare 
this  with  the  method  commonly  used  in  the  literature. 
Henceforth  referred  to  as  the  "incremental  distance"  method 
or  'method  1',  the  instantaneous  velocity  of  a  particle  is 
estimated  by  measuring  the  distance  a  particle  moves  between 
successive  motion  picture  frames,  and  dividing  this  distance 
by  the  period  T  between  successive  film  frames  (refer  to 
equation  (1.4)).  When  movement  in  the  y  coordinate  is  much 
smaller  than  that  in  the  x  coordinate,  as  it  is  in  the  case 
of  axonal  part icle  movement ,  equation  (1.4)  reduces  with 


little  error  to: 
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V((k+1)T)=  x( (k+1 ) I) -x(kT) 

T 

(3.1) 

This  is  equivalent  to  convolving  the  x  posi t ion- t ime  samples 
of  the  trajectory  with  a  two  term  approximation  to  an  ideal 
differentiating  filter.  The  magnitude  and  phase  responses 
can  be  shown  to  be  given  respectively  by: 


|H(f ) |=  (2/T)  sin  TT  fT 

if i  * 

1/2T 

(3.2) 

Z.H(f)=  IT  ( 1 -2f T) /2 

if i  * 

1/2T 

(3.3) 

This  frequency  response  is  compared  with 

that  of 

the  twenty 

one  term  filter  and  that  of 

an 

ideal 

low  pass 

differentiating  filter  in  figure  3.7. 

Simulated  trajectories  were  used  to  assess  the  error  of 
methods  1  and  2,  and  examine  the  effect  of  noise  in  the  data 
on  the  error  of  each. 
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PHASE  (DEGREES)  f]PGNITUDE  *10‘l 

0.0  22.5  45.0  67.5  90.0  0.0  15.0  30.0  45.0  60.0 
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(Hz) 


(a)  Magnitude  response  of  the  filters 


(b)  Phase  response  of  the  filters 


Figure  3.7  Comparison  of  the  filter  frequency  responses 
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Case  3.7(a) 

Three  simulated  trajectories  of  different  lengths  were 
analyzed  with  each  of  the  filters  to  estimate  the  errors 
inherent  in  each  method. 


Case  3.7(a)  results 


Method 

N 

(points) 

Vx  error 
(um/ sec) 

Vy  error 
(um/sec) 

Vd  error 
(um/sec) 

Mean 

S.D. 

Mean 

S.D. 

Mean 

S.D. 

1 

50 

-5.3E-3 

5.3E-2 

-4.8E-3 

3.9E-2 

-1 .  IE-2 

5.3E-2 

1 

100 

-6.2E-3 

4.5E-2 

-2.3E-3 

4. IE-2 

-8.7E-3 

4.5E-2 

1 

200 

-2.8E-3 

5.0E-2 

-2.5E-3 

4. IE-2 

-4.2E-2 

5.0E-2 

2 

50 

+5.7E-3 

1 .6E-2 

+1 .  6E-3 

8.4E-3 

+6.5E-3 

1 . 8E-2 

2 

100 

+7.7E-4 

1 .3E-2 

+9.5E-4 

1 .3E-2 

+1 . IE-3 

1 . 4E-2 

2 

200 

-2.4E-4 

1 .6E-2 

-1 .6E-4 

1 .3E-2 

-1 .2E-4 

1 . 6E-2 

Table  3.4 


Case  3.7(b) 

Gaussian  distributed  noise  of  mean=0.0  and  standard 
deviation  of  .13  urn  was  added  to  each  of  the  position 
coordinates  of  the  three  trajectories  previously  analysed. 
This  noise  level  was  equivalent  to  that  of  the  noisiest  real 
trajectory.  The  new  'noisy'  trajectories  were  analysed,  and 
the  errors  in  the  estimates  were  calculated. 


Case  3.7(b)  results 


Method 

N 

Vx  error 

Vy  error 

Vd  error 

(um/sec) 

(um/sec) 

(um/sec) 

(points) 

Mean 

S.D. 

Mean 

S.D. 

Mean 

S.D. 

1 

50 

-8.3E-3 

5.4E-1 

-3.3E-3 

5.2E-1 

-1 .3E-2 

5.3E-1 

1 

100 

-3.6E-3 

5.9E-1 

-5.3E-3 

5. IE-1 

-6.9E-3 

5.9E-1 

1 

200 

-2.6E-3 

5.7E-1 

-2.3E-3 

5.0E-1 

-3.7E-3 

5.7E-1 

2 

50 

+8.4E-3 

5. IE-2 

-2.3E-3 

5.3E-2 

+1 .OE-2 

5.8E-2 

2 

100 

+2.2E-3 

5.3E-2 

+8.7E-3 

5.9E-2 

+2.6E-3 

5.4E-2 

2 

200 

+6.9E-4 

4.5E-2 

+1 .3E-3 

5.2E-2 

+8.2E-4 

4.5E-2 

Table  3.5 
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VELOCITY  lU M/SEC. I  «ID‘  VELOCITY  lUft/SEC . I  ■10*' 

-10.0  0.0  10. 0  20.0  30.0  -LO.D  0.0  10. 0  20.0  30.0 
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0.0  10.0  20. □  30.0  40.0 

TIME  16ECJ 


TIME  tSEC] 


Actual  and  estimated  instantaneous  velocity 


Difference  between  actual  and  estimated  instantaneous  velocity 


Figure  3.8 


Estimation  of  Vx 


with  method  1 . 


a 


10.0 


ZD  .  D 

TIME  16ECJ 


Actual  and  estimated  instantaneous  velocity 


Difference  between  actual  and  estimated  instantaneous  velocity 


Figure  3.9  Estimation  of  Vx  with  method  2 
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Figure  3.10  Error  in  the  estimates  of  Vx 
as  a  function  of  noise  level. 
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Case  3.7(c) 

The  standard  deviation  of  the  Gaussian  noise  added  to 
the  three  trajectories  was  varied  to  determine  the 
relationship  between  noise  and  estimation  error  for  each  of 
the  methods.  The  method  error  for  estimates  of  Vx  is  given 
in  table  3.6  for  a  variety  of  cases. 


Case  3.7(c)  results 


N 

Noise 

(urn) 

method  error  (um/sec) | 

mean 

S.D. 

method  1 

method  2 

50 

0.0 

.065 

.31 

.036 

100 

0.0 

.065 

.28 

.033 

200 

0.0 

.065 

.26 

.032 

50 

0.0 

.130 

.53 

.052 

100 

0.0 

.130 

.55 

.056 

200 

0.0 

.130 

.54 

.049 

50 

0.0 

.260 

1.05 

.091 

100 

0.0 

.260 

1.14 

.087 

200 

0.0 

.260 

1.07 

.120 

Table  3.6 


It  should  be  noted  that  although  the  velocity  estimates 
and  associated  errors  discussed  throughout  this  work  are 
discrete  time  series,  for  reasons  of  clarity  they  are 
plotted  as  continuous  functions  of  time.  Figures  3.8  and  3.9 
illustrate  the  estimates  of  Vx  and  their  associated  errors 
obtained  with  the  two  methods  for  the  simulated  trajectory 
shown  in  figure  3.4.  The  effect  of  the  level  of  added  noise 
on  the  estimation  error  of  Vx  for  a  variety  of  cases  is 
tabulated  in  table  3.6.  Figure  3.10  illustrates  the  error 
with  different  noise  levels  for  N=100.  Method  2  was  found  to 
be  approximately  one  order  of  magnitude  less  sensitive  to 


noise  than  method  1. 
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3.8  Conclusions 

The  error  in  the  estimates  of  instantaneous  velocity 
obtained  with  the  method  developed  in  this  thesis  was 
approximately  one  order  of  magnitude  less  than  that  obtained 
with  the  incremental  distance  method  of  velocity  estimation. 
This  increase  in  accuracy  was  obtained  at  the  cost  of  the 
loss  of  ten  points  of  the  estimates  for  each  trajectory,  and 
an  increase  in  the  computing  time  required. 


CHAPTER  IV 
PROCEDURE 


As  reported  previously,  (Cooper  and  Smith,  1974),  16  mm 
motion  picture  films  at  framing  rates  of  3  and  4  frames/sec 
were  made  of  single  nerve  axons  of  Xenopus  1 aev i s  at  room 
temperature,  utilizing  darkfield  microscopy.  Two  populations 
of  moving  particles  were  observed  within  each  axon,  one 
consisting  of  particles  moving  somatopetal ly ,  and  one  of 
particles  moving  somatofugal ly . 

The  films  were  mounted  on  an  L-W  Photo  Optical  Data 
Analyzer  and  projected  onto  the  ground  glass  screen  of  a 
Photographic  Analysis  Ltd.  Digitizer.  Prior  to  filming  in 
all  cases  the  longitudinal  axis  of  the  axon  had  been 
oriented  such  that  on  the  film  records  it  was  parallel  to 
the  x  coordinate  axis  of  the  digitizer.  Segments  of  film 
were  repeatedly  viewed  to  locate  particle  trajectories 
satisfying  two  basic  requirements:  (1)  the  particle's 
position  could  be  clearly  discerned  for  a  large  number  of 
successive  frames  (preferably  >  100)  and  (2)  the  nerve 
preparation  did  not  move  appreciably  during  the  course  of 
the  trajectory.  Suitable  trajectories  were  recorded  on 
punched  paper  tape  by  making  frame  by  frame  readings  of  the 
selected  particle's  x  and  y  coordinates,  and  the  movie  frame 
number . 

The  trajectories  recorded  from  ten  normal,  healthy 
axons  (see  table  4.1)  were  selected  from  those  recorded  from 
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a  number  of  axons,  to  be  analyzed  in  detail.  The  optical 
quality  of  the  recorded  images  was  good,  and  there  was 
little  apparent  movement  in  each  case  of  the  axon  relative 
to  the  filming  apparatus.  One  hundred  thirty- three  of  the 
selected  trajectories  were  somatopetal,  thirty-nine  were 
somatofugal.  The  difference  in  population  sizes  reflected 
the  fact  that  a  greater  number  of  detectable  particles  move 
towards  the  cell  body  ( somatopetal ly)  than  away  from  it 
(somatofugal ly)  (Cooper  and  Smith, 1974). 


Tra iector ies  recorded 


Nerve 

Record  Des i gna  t i on 

Number  of  Trajectories  | 

Somatopetal 

Somatofugal 

A 

F,FF  Film  R1(a) 

10 

2 

B 

G,GG  Film  R1(b) 

10 

4 

C 

J,dd  Film  R48 

11 

2 

D 

M,MM  Film  R50 

10 

3 

E 

N,NN  Film  R51 

10 

7 

F 

0,00  Film  R52 

10 

1 

G 

P,PP  Film  R53 

10 

3 

H 

R,RR  Film  R35 

10 

3 

I 

S,SS  Film  R55 

42 

9 

J 

T , TT  Film  R54 

10 

5 

Total  Populations 

133 

39 

Table  4. 1 


A  typical  trajectory  is  shown  in  figure  4.1. 

The  spatial  character ist ics  of  particle  movement  were 
also  of  interest.  Axon  'I'  was  arbitrarily  selected  to  be 
analysed  in  particular  detail;  therefore,  more  trajectories 
(the  S  and  SS  series)  were  recorded  from  this  axon  than  from 
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Figure  4.1  Typical  particle  trajectory 
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Figure  A. 2  Axonal  structures  and  particle  trajectories. 
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the  others.  Also,  in  order  to  overlay  positional  data 
(i.e.  the  trajectories  and  the  related  velocity  components) 
it  was  necessary  to  show  that  the  positional  frame  of 
reference  (the  axon)  did  not  move  appreciably  relative  to 
the  digitizer.  Three  small,  separate,  apparently  stationary 
objects  were  located  in  this  axon  lying  within  the  field  of 
view  from  which  the  trajectories  had  been  recorded.  All 
three  objects  could  be  fairly  clearly  resolved  for  the 
entire  duration  of  the  film  segment.  As  it  was  assumed  that 
these  objects  were  indeed  stationary  relative  to  the  axon, 
their  positions  were  recorded  at  intervals  of  approximately 
300  frames  over  the  length  of  the  film  (approximately  2000 
frames)  to  obtain  an  estimate  of  the  amount  of  movement  of 
the  axon  relative  to  the  digitizer.  Apparent  movement  over 
an  1100  frame  segment  of  this  film  was  found  to  be 
particularly  small.  Over  this  segment  the  largest  ranges  in 
the  x  and  y  positions  of  the  three  objects  were  .52  urn  and 
.33  urn  respectively.  Therefore  the  maximum  apparent  movement 
of  the  axon  over  1100  frames  was  of  the  order  of  one 
particle  diameter,  and  was  considered  very  small  when 
compared  with  the  movement  of  particles  over  much  shorter 
film  segments.  Hence  it  was  concluded  that  spatial 
information  from  this  film  segment  could  be  "overlapped" 
with  little  positional  error  due  to  apparent  movement  of  the 
nerve.  In  addition,  position  readings  were  made  from  this 
film  segment  of  the  apparent  inner  surface  of  the  axon,  and 
of  long  stationary  rod-shaped  organelles  (probably 
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rod-shaped  mitochondria)  which  were  visible.  Figure  4.2  is  a 
plot  of  the  position  readings  of  the  inner  surfaces  of  the 
axonal  wall,  the  outlines  of  the  rod-shaped  organelles,  and 
the  particle  trajectories. 


< 


CHAPTER  V 
ANALYTICAL  METHODS 


5 . 1  Initial  process i na 

Despite  all  precautions,  the  position  samples 
occasionally  contained  gross  errors.  These  could  have  arisen 
from  dirt  on  the  digitizer  contact  wires,  or  errors  in  the 
digitizer  and  teletype.  Since  it  was  impossible  to  eliminate 
totally  this  type  of  error  at  the  source  ,  the  first 
procedure  in  the  analysis  of  the  data  attempted  to  identify 
and  correct  such  points.  It  was  essential  to  remove  these 
so-called  outliers  as  they  were  of  large  amplitude  and  could 
contain  a  substantial  portion  of  the  power  in  the  data. 

Chauvenet's  criterion  for  rejection  of  outliers  was 
used  to  identify  erroneous  readings  (Neville  and 
Kennedy , 1970) .  It  was  necessary  to  replace  such  points  with 
interpolated  values,  rather  than  merely  discarding  them, 
since  the  type  of  Fourier  analysis  used  assumed  the  data 
consisted  of  evenly  spaced  samples. 

The  x(kT)  and  y(kT)  sampled  position  sets  describing 
the  trajectories  were  stored  in  the  computer.  The  sets  were 
examined,  and  points  identified  as  outliers  were  discarded 
and  replaced  with  interpolated  values.  The  only  exceptions 
were  from  film  frames  in  which  the  particles  being  tracked 
had  been  observed  to  make  very  large,  sudden  movements.  Such 
movements  occurred  rarely,  but  were  due  to  "real"  movements 
of  the  particle,  not  instrument  or  observer  error,  and 
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therefore  were  not  replaced. 

A  third  sample  set;  the  radial  displacement  of  a 
particle  at  time  kT  from  its  initial  position  at  kT =0 , 
d(kT);  was  formed  from  the  corrected  x(kT)  and  y(kT)  data: 

d(kT)  =  /{x(kT)-x(0)}2+{y(kT)-y(0)}2 

k=0, 1 ,2. . .N-1 
N= number  of  samples 
T =sampl ing  period 

(5.1) 

The  time  series  x(kT),  y(kT),  and  d ( kT )  were  subsequently 
analyzed. 

5.2  Spectral  est imat ion 

Trend  removal  is  recommended  in  Fourier  analysis  to 
prevent  the  misinterpretation  of  spectrum  peaks  caused  by 
leakage,  and  to  prevent  leakage  from  masking  more  subtle 
phenomena.  Therefore,  linear  regression  lines  were 
calculated  for  the  sample  sets  x(kT),  y(kT),  and  d(kT);  and 
the  trends,  as  given  by  the  slopes  of  the  regression  lines 
(Vx,  Vy,  and  Vd  respectively)  and  their  axis  intercepts  (Bx, 
By,  and  Bd  respectively),  were  subtracted  from  each  of  them. 
The  results  were  trend -free  series  given  by: 
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X' (kT)=x(kT)-Vx  kT-Bx 

(5.2) 

y' (kT ) =y (kT) -Vy  kT-By 

(5.3) 

d/ (kT)=d(kT)-Vd  kT-Bd 

(5.4) 


k=0, 1 ,  .  . .  N-1 
N= number  of  samples 
T=  sampling  period 

Bx  =x  axis  intercept  of  linear  regression  line 

By  =y  axis  intercept  of  linear  regression  line 

Bd  =d  axis  intercept  of  linear  regression  line 


The  sets  x'(kT),  y#(kT),  and  d' (kT)  were  tapered  with 
the  Tukey  window  to  further  reduce  leakage,  producing 
tapered  trend- free  data  sets  x'MkT),  y'MkT),  and  d"(kT) 
respectively.  Computer  memory  space  available  with  the  power 
spectrum  estimation  program  limited  the  data  record  to  a 
maximum  of  512  points.  This  was  more  than  adequate,  allowing 
a  maximum  resolution  of  the  power  spectrum  of  about  .006  Hz 
(1/512T).  As  well,  none  of  the  data  sets  were  longer  than 
350  points,  well  below  this  limit.  Sufficient  zeroes  were 
added  to  the  end  of  the  data  sets  x'  '(kT),  y'  '(kT),  and 
d"  (kT)  to  pad  each  out  to  a  length  of  512  points.  This 
procedure  increased  the  resolution  of  the  resulting 
spectrum,  and  reduced  the  detrimental  effect  of  the  picket 
fence  phenomenon. 

The  windowed  trend-free  data  was  then  transformed  from 
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the  time  to  the  frequency  domain,  utilizing  the  FFT 
algorithm  as  implemented  by  Munro  (Munro, 1971 ) .  The 
transform  of  512  data  points  resulted  in  512  unique  Fourier 
coefficients,  which  sampled  the  frequency  space  from  zero  Hz 
up  to  one  half  the  sampling  frequency,  in  increments  of 
1/51 2T  Hz.  The  periodogram,  used  as  the  estimate  of  the 
power  spectrum,  was  calculated  from  the  Fourier  coefficients 
and  plotted. 

Using  these  procedures,  power  spectrum  estimations  were 
obtained  from  the  x77(kT),  y7/(kT),  andd77(kT)  data  sets 
obtained  from  one  hundred  and  seventy- two  trajectories 
tracked  in  a  total  of  ten  axons.  A  typical  power  spectrum 
consisted  of  several  peaks  below  .3  Hz.  The  power  above  .3 
Hz  was  very  small  compared  to  that  at  lower  frequencies.  It 


was  concluded  from  the  examination  of  a 

number 

of 

power 

spectra  that  the 

film  framing  rate 

(3 

Hz) 

had 

been 

sufficiently  high  to 

reduce  aliasing  error 

to 

a 

very 

smal  1 

value. 

The  power  at  low  frequency  in  y77(kT)  was  typically  one 
or  more  orders  orders  of  magnitude  smaller  than  that  in 
x77(kT),  and  the  d77(kT)  power  spectrum  was  almost  identical 
to  that  of  x77(kT).  As  the  x  component  of  particle  motion 
was  the  one  contributing  to  the  net  movement  of  the 
particles  along  the  axon,  it  was  the  power  spectra  of  the 
x77(kT)  series,  henceforth  simply  referred  to  as  the  7  power 
spectra7  ,  which  were  examined  in  detail. 

To  examine  the  differences,  if  any,  between  the  power 
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spectra  of  somatopetal  and  somatofugal  trajectories  the 
average  power  spectra  of  each  of  the  two  trajectory 
populations  was  calculated  from  the  power  spectra  of  the 
individual  trajector ies.  It  was  necessary  when  averaging  the 
power  spectra  to  introduce  a  scaling  factor  to  compensate 
for  the  varying  lengths  of  data  from  which  the  various 
spectra  had  been  obtained.  The  average  power  spectra  were 
calculated  with  the  following  equation: 


■§(mF) 


NS  A 
21  Kj  S j ( mF ) 
j  =  1 


(5.4) 


m=0, 1,2..  .N7- 

Ns -  the  number  of  spectra  averaged 
Kj-  the  scale  factor  N^/Nj 
N-p”  the  number  of  points  in  each  spectrum 
Nj -  the  number  of  samples  of  trajectory  j 


5.3  Determinat ion  of  f i 1  ter  cutoff  frequency 

Instantaneous  velocity  estimates  could  have  been 
obtained  by  wide  band  differentiation  of  the  sampled 
trajectory  data.  However,  with  this  method  the  effect  of  the 
very  weak  components  in  the  trajectory  at  higher  frequencies 
on  the  velocity  estimates  could  have  been  at  least  as 
significant  as  that  of  the  much  stronger  components  at  lower 
frequencies.  Before  proceeding  with  the  estimation  of 
particle  velocities  it  was  therefore  necessary  to  make  some 
assumptions  concerning  the  relative  significance  of  the 
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various  frequency  components  of  particle  motion. 

All  power  spectra  indicated  that  the  majority  of  power 
fell  in  the  frequency  band  from  0  to  .3  Hz.  At  frequencies 
above  .3  Hz  the  power  was  typically  three  orders  of 
magnitude  less  than  that  at  lower  frequencies.  These 
components  were  due  to  a  number  of  factors  including 
observer  error,  digitizer  error,  random  movements  of  the 
particle,  film,  microscope,  and  projector  movement.  Any 
component  in  the  data  due  to  measurement  error  was 
indistinguishable  from  that  due  to  "real"  random  movement  of 
the  particle. 

At  frequencies  below  about  .2  Hz  the  power  due  to 
measurement  error  was  insignificant  compared  to  the  power 
due  to  particle  motion,  and  hence  had  a  relatively 
insignificant  effect  on  the  estimates  of  instantaneous 
velocity.  However,  at  higher  frequencies  the  power  due  to 
measurement  error  was  estimated  to  be  at  least  as  great  as 
that  due  to  actual  particle  motion,  and  hence  could 
introduce  considerable  error  into  the  velocity  estimates. 
Consequently  it  was  necessary  to  remove  or  reduce  the 
effects  of  the  higher  frequency  components  on  the  estimates 
of  the  instantaneous  velocity.  A  low  pass  differentiating 
filter  was  used  in  the  estimation  procedure  to  accomplish 
this.  It  was  recognized  that  the  use  of  such  a  filter 
resulted  in  the  components  of  particle  motion  above  the 
cutoff  frequency  being  discarded  along  with  the  components 
due  to  measurement  error.  The  loss  of  these  components  of 
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particle  motion  was  unavoidable.  However,  their  loss  was 
considered  to  be  of  minor  consequence,  as  power  due  to  them 
was  so  much  less  than  the  power  at  frequencies  below  the 
filter  cutoff  frequency. 

If  the  average  power  level  of  the  higher  frequency 
components  had  been  consistent  from  trajectory  to 
trajectory,  the  filter  cutoff  frequency  for  each  trajectory 
analysis  could  have  been  specified  as  that  frequency  below 
which  the  power  in  the  spectrum  rose  above  a  level  based  on 
this  average  value.  Unfortunately,  the  high  frequency  power 
varied  with  each  trajectory,  as  one  could  track  a  small, 
sharply  defined  particle  which  remained  clearly  in  focus 
with  less  error  than  a  large,  fuzzy,  oval  particle  which 
dipped  in  and  out  of  focus.  Therefore,  the  filter  cutoff 
frequency  was  established  by  other  means. 

Ten  trajectories  were  digitized  and  their  power  spectra 
estimated.  The  percentages  of  total  power  falling  below  a 
number  of  cutoff  frequencies  were  calculated  for  a  group  of 
trajector ies ,  the  results  of  which  are  plotted  in  figure 
5.1.  As  the  cutoff  frequency  was  decreased,  the  percentage 
of  total  power  remaining  decreased  in  a  nonlinear  fashion. 
Decreasing  the  cutoff  frequency  from  .5  Hz  to  .2  Hz  did  not 
substantially  reduce  the  power  for  the  majority  of 
particles.  However,  as  the  cutoff  frequency  was  reduced 
below  .2  Hz  the  percentage  of  remaining  power  began  to  drop 
off  more  rapidly.  The  cutoff  frequency  was  conservatively 
chosen  to  be  .3  Hz,  and  for  the  sake  of  consistency  the  same 
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Figure  5.1  Effect  of  filter  cutoff  frequency  on  remaining 
signal  power. 
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cutoff  frequency  was  used  for  all  particles.  More  than  99% 
of  the  power  of  the  averaged  spectral  estimates  lay  below 
this  frequency. 

5.4  Velocity  est imat ion  by  low  pass  different iat ion 

The  optimum  procedure  for  instantaneous  velocity 
estimation  of  noisy  simulations  which  closely  resembled 
actual  trajectories  was  assumed  to  be  the  best  procedure  for 
velocity  estimation  of  actual  trajector ies.  Accordingly, 
estimations  of  the  average  and  instantaneous  velocities  of 
particles  were  carried  out  with  method  2,  developed  in 
Chapter  III. 

Estimates  of  the  x,  y,  and  d  components  of 
instantaneous  velocity  (henceforth  referred  to  as  Vx,  Vy, 
and  Vd  respectively)  were  computed  for  a  small  number  of 
trajectories.  In  addition,  the  magnitudes  and  phase  angles 
of  the  vector  sums  of  Vx  and  Vy  were  calculated.  The  various 
results  were  plotted  versus  time,  and  versus  position. 

It  was  essential  to  extract  the  information  of  most 
relevance  from  this  vast  quantity  of  data.  The  vector  sum  of 
Vx  and  Vy  was  found  to  be  of  little  value,  since  it  posed 
serious  problems  in  interpretat ion .  A  serious  problem 
associated  with  Vd  was  its  dependence  on  the  initial 
coordinates  of  each  set  of  trajectory  samples.  When 
examining  the  results  from  the  small  number  of  trajectories 
initially  analyzed,  it  was  seen  that  Vy  was  very  small  when 
compared  with  Vx,  and  in  most  cases  had  an  insignificant 
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contribution  to  the  net  movement  of  the  particle  (and 
therefore  presumably  an  insignificant  contribution  to  fast 
axoplasmic  transport).  As  Vx  lay  along  the  axis  along  which 
the  axonal  material  transport  occurs,  it  would  seem  to  be 
the  component  of  particle  motion  most  closely  related  to 
this  transport  process.  Therefore,  it  was  decided  that  Vx 
best  described  the  movement  of  particles  in  the  context  of 
the  physiological  process  of  interest  (axonal  transport),  in 
an  easily  interpretable  form.  It  was  Vx,  henceforth  referred 
to  as  the  "instantaneous  velocity",  which  was  therefore 
analyzed  in  detail. 

Each  particle  moved  predominantly  in  one  direction 
(i.e.  either  soma t ©petal ly  or  somatofugal ly ) ;  movement  and 
corresponding  velocity  in  the  predominating  direction  shall 
henceforth  be  referred  to  as  in  the  particle's  "major" 
direction.  Movement  and  corresponding  velocity  in  the 
direction  opposite  to  the  major  direction  shall  be  referred 
to  as  in  the  particle's  "minor"  direction.  The  instantaneous 
velocity  of  a  particle  (as  given  by  Vx)  moving  in  its  major 
direction  was  defined  to  be  positive;  and  in  its  minor 
direction,  negative. 

Instantaneous  velocity  estimations  were  computed  for 
one  hundred  and  seventy- two  trajectories  recorded  from  ten 
axons . 

5 . 5  Error  analysis 


5.5.1  Pos i t i on  record i nq  error 
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The  digitizer  error  was  found  to  correspond  to  ±  .015 
um  (±2SD,n=59)  at  the  magnification  used.  The  total  error 
in  estimating  the  position  of  a  moving  particle  was 
conservatively  estimated  to  be  ±  .08  um  (±  2SD, n=226) . 

5.5.2  Velocity  est imat ion  error 

The  error  in  the  velocity  estimation  method  was 
estimated  in  two  ways:  (a)  by  analyzing  simulations  whose 
character ist ics  closely  resembled  actual  trajector ies,  and 
(b)  by  redigitizing  the  film  records  of  previously  tracked 
trajectories  and  re-estimating  the  velocity.  To  determine 
the  consistency  of  the  estimation  method  the  estimates  thus 
obtained  were  compared  with  those  obtained  for  the  same 
particle  in  previous  trials. 

Analysis  of  simulated  particle  motion  (Chapter  III) 
indicated  that  the  x  and  y  components  of  particle 
instantaneous  velocity  could  be  estimated  with  a  maximum 
error  of  ±  .12  um/sec  (  db  2SD). 

The  difference  in  results  obtained  by  redigitizing  two 
trajectories  and  re-estimating  their  velocities  are 
tabulated  in  table  5.1. 
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Method 

N 

(points) 

Vx  difference 
(um/ sec) 

Vy  difference 
(um/sec) 

Vd  difference 
(um/sec) 

Mean 

S.D. 

Mean 

S.D. 

Mean 

S.D. 

1 

149 

-2.8E-3 

2. IE-1 

-4.2E-3 

1 .8E- 1 

-2.9E-3 

2.2E-1 

2 

149 

+2.5E-3 

3. IE-2 

+1 .3E-3 

2.9E-2 

+2.5E-3 

3. IE-2 

1 

88 

+7.8E-3 

2. IE-1 

-7.0E-3 

2.4E-1 

+7.2E-3 

2. IE-1 

2 

88 

+1 .4E-3 

1 . 9E-2 

+5.3E-4 

2.2E-2 

+1 .3E-3 

2. IE-2 

Table  5.1 


As  expected,  the  error  of  method  2  was  approximately  one 
order  of  magnitude  less  than  that  of  7  method  1' ,  and  was 
well  within  the  limits  predicted  by  simulation  analysis. 
Estimates  of  the  instantaneous  velocity  (given  by  Vx,  as 
previously  discussed)  obtained  by  digitizing  a  trajectory 
twice  and  repeating  the  velocity  estimation  using  the  two 
different  methods  may  be  seen  in  figure  5.2.  Also  shown  is 
the  difference  between  the  estimates  for  each  of  the 
methods. 

5.5.3  Effect  of  ignoring  y  and  z  components  of  pos i t i on 

Throughout  this  work  the  movement  of  particles  in  the 
axon  was  treated  as  being  two-dimensional.  Movement  in  the  z 
direction;  that  is,  in  the  direction  perpendicular  to  the 
microscope's  plane  of  focus,  was  entirely  ignored,  since  it 
was  impossible  to  measure. 

It  was,  however,  possible  to  estimate  the  magnitude  of 
this  movement  through  two  approaches:  (a)  consideration  of 
the  optical  character ist ics  of  the  cinemicrophotographic 
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(c)  Difference  between  estimates  with  Method  1  (d)  Difference  between  estimates  with  Method 
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apparatus,  and  (b)  consideration  of  the  magnitude  of  the 
movement  in  the  y  coordinate. 

It  has  been  estimated  (Smith,  personal  communication) 
that  the  depth  of  focus  of  the  optical  system  used  in  this 
work  through  which  these  axonal  particles  can  be  discerned 
was  approximately  1.5  urn.  That  is,  for  a  particle  to  be 
resolved  for  the  entire  length  of  its  trajectory,  its  total 
movement  in  the  z  direction  must  be  less  than  1.5  urn.  This 
movement  was  small  when  compared  with  the  typical  trajectory 
length  of  15  to  30  urn  and  implied  that  the  average  z 
velocity  component  of  such  particles  was  very  small. 

To  obtain  a  more  quantitative  estimate  of  movement  in 
the  z  direction,  movements  of  particles  in  the  y  coordinate 
were  considered.  It  has  been  seen  in  this  work  that  the 
displacements  (and  corresponding  instantaneous  velocity 
components)  in  the  y  coordinate  were  almost  always  at  least 
one  order  of  magnitude  smaller  than  those  in  the  x 
coordinate.  Since  the  y  coordinate  arbitrarily  determined 
the  viewing  plane  through  the  axon  it  was  reasonable  to 
assume  that  movement  of  particles  in  the  z  direction,  also 
perpendicular  to  the  longitudinal  axis  of  the  axon,  was  of 
the  same  order  of  magnitude. 

In  this  work  the  x  component  of  velocity  was  used  as 
the  "particle  instantaneous  velocity".  The  error  due  to 
ignoring  the  y  and  z  components  of  velocity  was  considered 
very  small,  and  was  conservatively  estimated  to  be 
approximately  ±  .05  um/sec  (  di  2SD). 
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5.5.4  Conclusions  of  error  analysis 

The  error  in  instantaneous  velocity  estimates  as 
determined  by  retracking  and  reanalyzing  several 
trajectories  was  approximately  one  half  that  of  the  error  as 
predicted  by  the  analysis  of  simulated  tra jector ies .  If  the 
worst  case  is  taken,  and  the  error  due  to  ignoring  the  y  and 
z  components  of  velocity  is  included,  the  error  in  the 
particle  instantaneous  velocity  estimates  used  in  this  work 
is  ±  .17  um/sec  (  ±,  2SD). 

5.6  Spatial  Analysis 

Many  particle  trajectories  could  be  recorded  from 
observation  of  an  axon  at  a  single  location.  In  addition, 
the  use  of  improved  techniques  in  nerve  preparation, 
cinemicrophotography ,  and  analytical  techniques  enabled 
accurate  estimation  of  the  instantaneous  velocities  of  a 
large  number  of  particles  in  a  single  axon.  These 
improvements  in  the  quality  and  quantity  of  data  available 
facilitated  the  exam inat ion  of  the  spatial  relationships 
between  particle  velocities.  The  scope  of  such  spatial 
analysis  is  vast  indeed,  and  the  topic  was  therefore  given 
only  preliminary  examination.  It  is  hoped,  however,  that 
this  work  will  suggest  approaches  for  future  trajectory 
analysis. 

The  topics  of  mapping  and  spatial  analysis  have  been 
widely  examined  in  the  geographical  and  geophysical 
literature.  However,  many  of  the  principles  and  techniques 
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thus  developed  have  a  wider  scope  and  are  suitable  for  use 
in  other  disciplines. 

Associated  with  each  instantaneous  velocity  estimate 
were  corresponding  x  and  y  position  readings.  The  data  could 
be  considered  to  be  three  dimensional,  and  was  similar  in 
concept  to  the  data  used  in  producing  topographical  contour 
maps.  The  x  and  y  position  readings  corresponded  to  the 
geographical  coordinates  of  the  survey  points;  and  the 
instantaneous  velocity  estimates  to  the  elevation  estimates. 
SURFACE  II  ,  a  computer  software  system  initially  developed 
for  use  in  topographical  mapping,  was  found  to  be  suitable 
for  use  in  examining  the  velocity-position  data. 

SURFACE  II  was  developed  by  the  Kansas  Geological 
Survey  in  the  early  1970's.  The  system  is  quite  general,  and 
will  display  the  form  of  any  variable  characterised  by 
values  at  coordinates  defined  by  two  other  variables.  The 
only  inherent  restrictions  are  that  the  coordinate  variables 
must  be  orthogonal,  and  the  mapped  variable  must  be  single 
valued.  A  basic  form  of  graphic  display  produced  by  SURFACE 
II  is  a  contour  map,  a  plot  of  the  two  coordinates  on  which 
values  of  the  third  variable  are  defined  by  lines  of  equal 
value. 

The  basis  for  all  operations  in  SURFACE  II  is  a 
rectangular  grid  of  values,  which  is  a  numerical 
representat ion  of  the  surface  to  be  displayed.  The  grid  is  a 
numerical  matrix  whose  elements  are  measurements  of  the  z 
variable  made  at  points  identified  by  the  two  'geographic7 
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coordinates,  x  and  y.  The  contour  plot  representat ion  of 
this  matrix  is  obtained  by  linearly  interpolating  between 
grid  nodes  to  locate  the  points  where  a  contour  line  of 
specified  value  will  cross  the  edge  of  a  grid  cell.  The 
string  of  successive  x  and  y  coordinates  of  these 
intersections  defines  the  contour  line  that  will  be  drawn. 

The  velocity-position  data  consisted  of  points 
scattered  irregularly  over  the  'map7  area.  Since  the 
graphics  display  required  a  regular  grid,  the  first 
procedure  was  to  estimate  such  a  grid  from  the  irregularly 
spaced  data. 

There  are  two  basic  classes  of  methods  for  such  grid 
estimation:  global  fit  methods,  and  local  fit  methods.  With 
global  fit  methods  all  data  points  are  used  in  the 
estimation  of  each  grid  node  value,  whereas  local  fit 
methods  assume  that  only  the  points  in  the  "neighbourhood" 
of  each  node  should  be  used  in  the  estimation  of  its  value. 
Global  fit  methods  approximate  the  data  values  with  three 
dimensional  functions  and  use  interpolation  to  estimate  the 
values  at  the  grid  nodes;  local  fit  methods  estimate  the 
node  values  with  distance-weighted  averages  of  the 
"nearest -neighbour"  data  points  of  each  node.  Local  fit 
methods  are  generally  considered  to  be  the  most  appropriate 
for  estimating  points  on  complex  surfaces  (Sampson, 1975) , 
and  were  judged  the  most  appropriate  for  use  in  the 
examination  of  spatial  character ist ics  of  particle 


veloci t ies. 
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The  "nearest -neighbours"  of  a  grid  node  may  be  defined 
by  a  number  of  search  techniques,  some  of  which  may  provide 
"better"  results  than  others  in  certain  situations 
(primarily  reflecting  the  nature  of  the  data  point 
distribution).  For  example,  some  methods  may  not  be  suitable 
when  data  points  are  collected  along  lines  or  traverses,  as 
in  the  case  of  trajectory  analysis.  In  such  a  case,  a  simple 
nearest  neighbour  search  may  find  the  specified  number  of 
points  for  an  estimation  entirely  from  a  single  trajectory. 
Such  an  estimate  would  be  constrained  in  only  one  direction. 
The  quadrant  search  technique,  among  others,  is  recommended 
in  such  cases  to  ensure  some  radial  control  in  the 
estimation  procedure  (Sampson,  1975). 

A  quadrant  search  procedure  divides  the  area  around  a 
grid  node  into  four  equal  segments.  The  nearest  n  data 
points,  where  n  is  one  of  the  search  procedure  parameters, 
are  found  in  each  segment,  within  a  maximum  search  radius 
specified  by  another  of  the  search  parameters.  If  the 
distance  to  the  nearest  data  point  exceeds  a  specified 
distance,  the  search  fails  and  the  grid  node  is  not 
estimated.  Should  less  than  a  specified  number  of  quadrants 
contain  at  least  one  data  point  within  the  maximum  search 
radius,  the  search  will  fail  and  the  grid  node  will  not  be 
estimated.  Assuming  the  search  has  been  successful,  the 
value  at  the  grid  node  is  estimated  by  calculating  a 
distance-weighted  average  of  the  data  points  found. 

The  x-y  coordinate  positions  and  corresponding 
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instantaneous  velocities  of  forty  two  somatopetal ly  moving 
particles,  recorded  in  axon  I,  were  stored  on  paper  tape. 
This  data  was  input  to  the  University  of  Alberta  Computing 
Services  Amdahl  470  computer,  and  with  the  use  of  the 
SURFACE  II  software  package,  a  number  of  contour  maps  were 
produced,  using  different  algorithms  and  nearest -neighbours 
search  techniques  to  estimate  grids  from  the 
velocity-position  data. 

Rectangular  grid  matrices  were  estimated  from  the  4375 
irregularly  spaced  x-y  coordinate  positions  and 
corresponding  instantaneous  velocities.  The  matrices 
encompassed  an  area  of  18  urn  by  38  urn,  and  were  evenly 
divided  into  684  cells,  each  1  urn  by  1  urn.  Selection  of  the 
cell  size  involved  a  tradeoff  between  minimizing  computing 
time  (i.e.  using  a  small  number  of  large  cells),  and 
maximizing  spatial  resolution  (i.e.  using  a  large  number  of 
sma 11  cells). 

The  quadrant  search  technique  was  used  to  estimate  the 
values  at  the  grid  nodes  from  the  velocity-position  data. 
Determining  suitable  values  for  the  search  parameters 
required  considerable  planning,  and  a  certain  amount  of 
trial  and  error.  The  estimation  procedure  in  effect  low  pass 
filtered  the  data.  Therefore,  the  area  over  which  the  data 
points  were  averaged  had  to  be  large  enough  to  include 
samples  from  several  tra jector ies ,  and  yet  small  enough  to 
be  able  to  resolve  localized  variations  in  velocity. 

The  values  used  in  the  final  analysis  specified  a 
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maximum  search  radius  of  1.5  urn.  At  least  one  point  had  to 
lie  within  1.0  urn  of  the  grid  node.  A  minimum  of  3,  and  a 
maximum  of  12,  points  were  used  to  estimate  the  value  at 
each  node.  The  points  selected  by  the  search  algorithm  were 
weighted  by  the  function  1/D  (where  D  was  the  distance  from 
the  point  to  the  grid  node),  and  the  value  at  the  node 
obtained  by  averaging  the  weighted  values. 


CHAPTER  VI 


RESULTS  OF  THE  ANALYSIS  OF  THE  MOTION 
OF  A  GROUP  OF  INTRA-AXONAL  PARTICLES 

6 . 1  Spectral  est imat ion 

Examination  of  a  number  of  power  spectra  indicated 
appreciable  power  at  low  frequencies.  Above  approximately  .1 
Hz  the  power  in  all  cases  rapidly  decreased  by  more  than  two 
orders  of  magnitude  to  an  insignificant  level.  At  higher 
frequencies  the  power  presumably  arose  from  random 
components  in  the  data,  including  all  measurement  error,  and 
any  "real"  random  movements  of  the  particles  which  were 
indistinguishable  from  random  error. 

Two  averaged  power  spectra  were  calculated  from  the 
power  spectra  obtained  from  the  two  trajectory  populations 
(i.e.  133  somatopetal  trajector ies,  and  39  somatofugal 
trajector ies) ,  and  are  plotted  in  figure  6.1.  The  shapes  of 
the  averaged  power  spectra  for  the  two  populations  were 
similar.  At  low  frequencies  the  average  power  of  the 
somatopetal  trajectories  was  somewhat  more  than  twice  that 
of  the  somatofugal  trajector ies,  although  at  higher 
frequencies  the  average  power  of  the  two  populations  was  the 
same. 

6 . 2  Average  veloci t ies 

The  x,  y,  and  d  component  average  velocity  estimates 
(Vx,  Vy”,  and  Vd  respectively)  were  given  by  the  slopes  of 
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(a)  Averaged  power  spectrum  obtained  from  133  somatopetal  trajectories. 


(b)  Averaged  power  spectrum  obtained  from  39  somatofugal  trajectories. 


Figure  6.1  Averaged  power  spectra 
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the  linear  regression  lines  of  sampled  sets  x(kT),  y(kT), 
and  d(kT)  respectively  on  time.  The  results  for  172  particle 
trajectories  recorded  from  a  total  of  ten  axons  are 
tabulated  in  tables  6.1  and  6.2.  Histograms  of  the  average 
velocities  for  the  somatofugal  and  somatopetal  populations 
may  be  seen  in  figure  6.2. 
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Average  velocity  statist ics  by  axon 


Axon 

Somatopetal 

Trajector ies 

Somatofugal 

Trajector ies 

n 

Vx  Mean 

S.D. 

Range 

n 

Vx  Mean 

S.D. 

Range 

(um/sec) 

min 

max 

(um/sec) 

min 

max 

A 

10 

.49 

.38 

.06 

1.37 

2 

.98 

.33 

.75 

1.21 

B 

10 

.93 

.41 

.28 

1.42 

4 

1 . 10 

.62 

.58 

1 .89 

C 

11 

.83 

.41 

.30 

1.44 

2 

.67 

.06 

.63 

.71 

D 

10 

1.07 

.40 

.57 

1 .91 

3 

.56 

.  14 

.41 

.68 

E 

10 

.86 

.34 

.35 

1.43 

7 

.67 

.06 

.59 

.75 

F 

10 

1.12 

.35 

.71 

1.66 

1 

1.14 

G 

10 

1.16 

.43 

.39 

1 .68 

3 

.89 

.42 

.56 

1.36 

H 

10 

.91 

.32 

.45 

1 .41 

3 

1.01 

.60 

.63 

1.71 

I 

42 

.81 

.32 

.23 

1 .59 

9 

1.04 

.60 

.27 

1 .93 

J 

10 

.98 

.42 

.32 

1 .81 

5 

1.07 

.72 

.25 

1 .95 

Table  6.1 


Pooled  average  velocity  statist ics 


Component 

Somatopetal 

Somatofugal 

Mean 

S.D. 

Range 

Mean 

S.D. 

Range 

(um/sec) 

min 

max 

(um/sec) 

min 

max 

Vx 

.91 

00 

CO 

■ 

.06 

1 .91 

.93 

.47 

.25 

1 .95 

V  y 

.00 

.11 

-.29 

.47 

.02 

.10 

-.27 

.32 

Vd 

.91 

.38 

.06 

1.96 

.93 

.47 

.25 

1  .94 

n=  133 

n=39 

Table  6.2 
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The  somatofugal  and  somatopetal  populations  of 
trajectory  average  velocity  were  examined  to  determine  the 
differences,  if  any,  between  them.  Using  the  F  test,  the 
variances  of  the  somatofugal  and  somatopetal  populations 
were  found  to  be  significantly  different  at  the  5%  level. 
Therefore,  to  examine  the  difference  in  the  population 
means,  it  was  necessary  to  use  a  test  which  does  not  assume 
that  the  two  population  variances  are  the  same.  Using  such  a 
test  (p. 142,Nevi lie  and  Kennedy , 1970) ,  the  difference  in 
mean  velocities  between  the  two  populations  was  found  not  to 
be  statistically  significant  at  the  5%  confidence  level. 

6 . 3  Instantaneous  veloci t ies 

Examination  of  plots  of  particle  instantaneous  velocity 
versus  time  showed  instantaneous  velocities  which  smoothly 
varied  (fig.  6. 3(a), (b)).  Intuitively  these  results  seemed 
suspect,  since  the  jerky  or  saltatory  movement  of  axonal 
particles  is  a  widely  described  phenomenon.  If,  however, 
plots  of  instantaneous  velocity  versus  distance  moved  were 
examined,  saltatory  movements  were  clearly  evident  (fig. 
6. 3(c), (d)).  Thus  it  was  concluded  that  although  the 
particles  move  in  a  saltatory  manner  with  respect  to 
position,  their  velocities  are  smooth,  continuously  varying 
functions  with  respect  to  time. 

To  examine  the  distribution  of  the  instantaneous 
velocity  estimates  for  each  of  the  two  particle  populations 
the  composite  frequency  distribution  of  instantaneous 
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Figure  6.3 
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NEGATIVE  -  «L  S  X  POSITIVE  -  JJ.  5  X. 
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Figure  6 . 4 
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(a)  Somatopetal  trajectories 


(b)  Somatofugal  trajectories 


Figure  6.5  Distribution  of  the  standard  deviation  of  the 
instantaneous  velocity. 
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Figure  6.6  Standard  deviation  of  the  instantaneous  velocity 
plotted  versus  the  mean  velocity. 
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velocity  for  each  was  computed  (figure  6.4).  It  was 
necessary  when  computing  the  composite  distributions  to 
normalize  the  contributions  of  each  trajectory,  since  in 
general  slower  moving  particles  contributed  more 
instantaneous  velocity  estimates.  This  was  done  by  weighting 
the  contributions  of  each  trajectory  to  each  "bin"  of  the 
composite  histogram  by  1/N,  where  N  is  the  total  number  of 
estimates  contributed  by  the  trajectory. 

The  distribution  of  the  standard  deviation  of 
instantaneous  velocity  of  each  particle  about  its  respective 
mean  value  for  the  two  populations  may  be  seen  in  figure 
6.5.  The  mean  of  the  somatopetal  distribution  is 
substantially  greater  than  that  of  the  somatofugal 
distribution,  implying  that  on  the  average  instantaneous 
velocity  of  somatopetally  moving  particles  varies  more 
widely  than  that  of  somatofugally  moving  particles. 

To  examine  the  relationship,  if  any,  of  the  amount  by 
which  the  velocity  of  a  particle  varied  about  its  mean  value 
to  that  mean  value,  the  standard  deviations  of  instantaneous 
velocity  for  the  172  trajectories  were  plotted  versus  their 
mean  values  (fig  6.6).  The  line  of  linear  regression  of  the 
standard  deviations  on  the  mean  velocities  was  computed,  as 
was  the  correlation  coefficient.  These  results  indicated  no 
correlation  between  the  amount  by  which  the  instantaneous 
velocity  of  particles  varied  about  their  mean  values  and 
those  mean  values. 

In  the  search  for  patterns  or  groupings  in  the  movement 
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of  particles  it  was  essential  to  quantify  their  motion  and 
to  extract  from  the  vast  quantity  of  instantaneous  velocity 
data  available  a  small  number  of  descriptive  parameters. 

An  index  of  the  net  progress  of  particles,  known  as  the 
vector-scalar  ratio  (VSR),  has  been  defined  by  Breuer 
(Breuer  et  a  1., 1975).  The  VSR  is  the  ratio  of  the  net 
distance  a  particle  moves  divided  by  its  total  path  length. 
The  VSR  may  be  considered  as  a  measure  of  "efficiency"  of 
particle  movement.  For  example,  a  particle  which  never 
reversed  its  direction  of  movement  would  have  a  VSR  of  1;  at 
the  other  extreme,  a  particle  which  made  no  net  movement 
would  have  a  VSR  of  0. 

The  VSR  was  calculated  for  each  of  the  133  somatopetal 
and  39  somatofugal  trajectories.  Histograms  of  the  VSR's  are 
plotted  in  figure  6.7.  In  figure  6.8  the  trajectory  VSRs  are 
plotted  versus  the  trajectory  mean  velocities,  for  the  two 
populations.  The  line  of  linear  regression  of  the  trajectory 
VSR's  on  the  mean  velocities  was  computed,  as  was  the 
correlation  coefficient.  Correlation  was  significant  at  the 
99%  level,  and  it  was  concluded  that  the  smaller  the  average 
velocity  of  a  given  particle,  the  smaller  its  VSR. 

Another  parameter  of  interest  was  the  percentage  of 
time  a  particle  spent  moving  in  the  minor  direction 
(i.e.  the  percentage  of  velocity  estimates  which  were 
negative  for  a  given  trajectory).  Figure  6.9  illustrates  the 
distribution  of  this  parameter  for  the  two  trajectory 
populations.  In  figure  6.10  the  percentages  of  time 
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Distribution  of  VSR  for  somatopetal  trajectories. 


(b)  Distribution  of  VSR  for  somatofugal  trajectories. 
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Figure  6.8  Relationship  of  trajectory  VSR  to  average  velocity. 
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(a)  Distribution  of  percent  of  time  spent  moving  in  minor 
direction  for  somatopetal  trajectories. 
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(b)  Distribution  of  percent  of  time  spent  moving  in  minor 
direction  for  somatofugal  trajectories. 
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Figure  6.10  Relationship  of  time  spent  travelling  in  minor 
direction  to  average  velocity. 
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Figure  6.11  Relationship  of  trajectory  VSR  to  the  amount  of  time  spent 
travelling  in  the  minor  direction 
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particles  spent  moving  in  the  minor  direction  were  plotted 
versus  their  mean  velocities.  The  line  of  linear  regression 
of  the  percentage  of  time  in  the  minor  direction  on  the  mean 
velocity  was  computed,  as  was  the  correlation  coefficient. 
The  correlation  between  the  variables  was  significant  at  the 
99%  level,  and  it  was  concluded  that  the  smaller  the  mean 
velocity  of  a  given  particle,  the  greater  the  percentage  of 
time  it  was  likely  to  spend  moving  in  the  minor  direction. 

Figure  6.11  illustrates  the  relationship  between  the 
net  movement  of  particles  (as  measured  by  the  VSR),  and  the 
amount  of  time  the  particles  spent  moving  in  the  minor 
direct  ion . 

6.4  Spatial  aspects  of  part icle  mot  ion 

Figures  6.12(a)-(d)  are  plots  of  particle  trajectories 
and  their  corresponding  velocities  for  four  pairs  of 
trajectories.  Each  pair  represents  two  particles  whose  paths 
overlapped  in  space  for  a  portion  of  their  trajectories.  In 
each  pair  of  observations  a  similarity  was  apparent  in  the 
velocities  over  those  portions  of  the  trajectories  which 
overlapped.  Although  similarities  in  velocities  of 
overlapping  trajectories  were  not  apparent  in  all  cases, 
they  occurred  with  sufficient  frequency  to  warrant  further 
invest igat ion. 

It  can  be  hypothesized  that  particles  in  the  same 
spatial  location  in  the  axon  move  with  similar  velocities. 
This  would  be  expected  if  the  force  generating  mechanism 
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Figure  6.12  Overlapping  trajectory  pairs  and  their  corresponding 
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Figure  6.13  Contour  map  representing  spatial  distribution  of  particle  position 
measurements  and  corresponding  instantaneous  velocity  estimates. 
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Figure  6.14  Contour  map  representing  the  spatial  distribution  of  particle  instantaneous  velocities. 
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responsible  for  particle  motion  consisted  of  discrete, 
"stationary"  sites  distributed  throughout  the  axon.  It  would 
also  be  the  case  if  the  particle  axonal  velocities  were 
significantly  related  to  local  variations  in  the  properties 
of  the  axoplasm.  This  hypothesis  was  investigated  with  the 
use  of  contour  plots. 

Figures  6.13  and  6.14  were  obtained  by  analyzing  4375 
estimates  of  velocity  obtained  from  42  somatopeta 1 ly  moving 
particles  in  axon  'I'.  Figure  6.13  is  a  graphical 
representation  of  the  distribution  of  data  points,  i.e  the 
geographical  density  of  the  velocity  estimates.  It  was 
obtained  by  computing  the  number  of  velocity  estimates  lying 
within  .707  urn  of  each  grid  node,  and  then  computing  and 
plotting  a  contour  map  of  the  number  of  points  around  each 
grid  node  versus  x  and  y  position.  The  search  radius  for  the 
density  estimates  was  chosen  to  be  .707  urn  to  ensure  that 
all  the  area  between  the  grid  nodes  (spaced  1  urn  apart)  was 
covered,  with  minimal  overlap. 

The  velocity-position  contour  plot  shown  in  figure  6.14 
attempts  to  display  graphically  the  combined  velocity 
results  for  forty  two  somatopetally  moving  particles.  Should 
the  velocity  of  each  particle  be  totally  unrelated  to  that 
of  other  particles  moving  in  the  same  region  in  the  axon, 
the  velocity  estimates  at  the  grid  nodes  should  be  quite 
similar.  The  contour  plot  would  therefore  consist  of  only  a 
few  contour  lines.  This  is  not  the  case  however,  implying  a 
correlation  between  position  and  velocity. 
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As  can  be  seen  from  figure  6.14,  the  velocities  of 
particles  averaged  in  various  regions  of  the  axon,  ranged 
from  0.0  um/sec  to  2.0  um/sec.  Areas  of  similar  velocities 
varied  considerably  in  size,  and  were  longitudinally 
arranged  with  respect  to  the  axon.  There  were  no  obvious 
relationships  between  these  areas  and  optically  discernible 
axonal  structures  such  as  the  rod  shaped  mitochondria. 

A  similar  examination  of  spatial  aspects  of  somatofugal 
particle  motion  could  not  be  conducted  since  insufficient 
measurements  were  available  in  individual  axons. 


CHAPTER  VII 


DISCUSSION  AND  CONCLUSIONS 


7 . 1  Introduct ion 

The  preceding  chapters  have  outlined  the  analysis  of 
axonal  particle  movement  by  a  variety  of  methods.  This 
chapter  shall  compare  the  results  thus  obtained  with  those 
of  other  workers  in  the  field,  and  shall  attempt  to  draw 
some  conclusions  on  the  nature  of  axonal  particle  motion. 

7 . 2  Spectral  analysis 

In  their  analysis  of  axonal  particle  movement,  Forman's 
group  (Forman  et  a  1,1977a)  searched  for  periodicity  in  the 
instantaneous  velocity  estimates  of  61  individual  particles 
by  means  of  autocorrelat ion  and  spectral  analysis.  They 
reported  that  the  results  were  inconclusive,  with  no 
conspicuous  or  consistent  patterns  emerging.  In  retrospect 
this  is  not  suprising,  since  the  instantaneous  velocity 
estimates  had  been  obtained  with  the  incremental  distance 
method,  and  hence  contained  considerably  amplified  noise. 
The  power  at  higher  frequencies  due  to  the  noise  in  the 
velocity  estimates  could  be  at  least  as  great  as  that  at 
lower  frequencies  due  to  both  noise  and  actual  particle 
motion.  Any  underlying  periodic  components  at  low 
frequencies  would  not  be  obvious. 

In  this  work  the  power  spectra  of  172  particle 
trajectories  were  analyzed.  Any  consistent  periodic 
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component  in  the  motion  of  the  particles  would  be  evidenced 
by  peaks  in  the  autopower  spectra  of  the  trajectory  and 
velocity  estimates.  The  trajectory  power  spectra  were 
estimated  in  the  manner  previously  described. 

The  considerable  differences  in  the  spectral  estimates 
obtained  with  the  two  different  methods  (i.e.  Forman's  vs 
McLeod's)  may  be  seen  for  two  trajectories  in  figures 
7. 1 (a)-(d) .  Figures  7.1(a)  and  (b)  are  spectral  estimates 
obtained  by  analyzing  two  trajector ies.  Figures  7.1(c)  and 
(d)  were  obtained  by  estimating  the  instantaneous  velocities 
of  the  same  two  trajectories  using  the  incremental  distance 
method,  and  then  analyzing  the  velocity  estimates.  The  power 
at  higher  frequencies  due  to  the  accentuated  noise  in  the 
velocity  estimates  is  evident. 

The  trajectory  power  spectra  suggested  a  general 
pattern  in  the  motion  of  particles.  Particle  motion  was 
consistently  composed  of  low  frequency  components,  typically 
in  the  range  of  from  0  to  .2  Hz.  Above  .2  Hz  the  magnitude 
of  the  frequency  components  was  much  smaller.  Thus  it  was 
concluded  that  the  particles  moved  with  an  underlying  low 
frequency  periodic  motion,  whose  period,  although  varying 
from  particle  to  particle,  was  consistently  longer  than  5 
seconds.  It  was  also  evident  that  a  filming  rate  of  3  fps 
was  more  than  sufficient  to  record  the  major  details  of 
part icle  mot  ion . 

The  peak  power  of  the  averaged  power  spectra  of  the 
somatopetal  trajectories  was  found  to  be  somewhat  more  than 
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(c)  Trajectory  1  -  Velocity  power  spectrum  (d)  Trajectory  2  -  Velocity  power  spectrum 

as  obtained  with  Forman's  procedure.  as  obtained  with  Forman's  procedure. 
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twice  that  of  the  somatofugal  trajector ies.  This  implies 
that  the  instantaneous  velocity  of  a  somatopetal  particle 
tends  to  vary  more  widely  about  its  average  value  than  does 
that  of  a  somatofugally  moving  particle.  Examination  of  the 
distribution  of  the  standard  deviations  of  the  instantaneous 
velocity  for  the  two  populations  (fig  6.5)  leads  to  a 
similar  conclusion. 

The  power  spectra  of  the  x  component  of  position  was 
used  as  the  "trajectory  power  spectra".  Since  the  power  in 
the  y  and  presumably  the  z  components  was  two  or  more  orders 
of  magnitude  less  than  that  in  the  x  component  in  the  region 
of  interest,  the  consequences  of  ignoring  them  were  not 
significant . 

7.3  Particle  average  veloci t ies 

The  average  velocities  of  somatofugal  and  somatopetal 
particles  were  similar  to  the  values  reported  by  other 
workers  with  similar  preparations  (Cooper  and  Smith, 1974; 
Smith  and  Koles,1976;  Forman  et .  al.,1977a,b)  and  observed 
in  intact  myelinated  axons  (Hammond, 1977a) .  Although  the 
means  of  the  two  populations  of  particle  average  velocities 
were  not  significantly  different,  the  variance  of  the 
somatofugal  population  was  significantly  greater  than  that 
of  the  somatopetal  population.  This  was  consistent  with  the 
results  reported  by  Hammond  (Hammond, 1977a ) . 


7 . 4  Velocity  est imat ion  method 
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Analyses  of  both  simulated  and  actual  particle 
trajectories  have  clearly  shown  the  advantages  of  the  method 
developed  in  this  thesis  for  the  estimation  of  particle 
instantaneous  velocity  over  the  method  found  in  the 
literature.  It  was  felt  that  the  reduction  of  the  estimation 
error  by  at  least  one  order  of  magnitude  more  than 
compensated  for  the  loss  with  this  method  of  approximately 
10%  of  the  results.  With  more  accurate  estimation  of 
instantaneous  velocity  can  come  more  sophisticated  and  more 
comprehensive  analysis  of  particle  instantaneous  velocity. 

As  discussed  previously  for  practical  reasons  in  this 
work  the  x  component  of  velocity  was  used  as  the  particle 
instantaneous  velocity,  and  the  y  and  z  components  were 
ignored.  The  error  due  to  ignoring  these  components  was 
found  to  be  very  small.  This  approach  must  be  used  with  care 
however,  since  in  a  very  small  number  of  cases  (2-4%) 
particles  were  observed  to  make  substantial  movements  in  the 
y  direction.  Workers  in  the  future  may  wish  to  examine  such 
trajectories  as  special  cases  using  more  detailed  analysis. 
Substantial  movement  in  the  z  direction  was  not  a  problem, 
as  in  such  cases  the  particle  would  leave  the  thin  focal 
plane  of  the  filming  apparatus  (disappearing  from  view)  and 
hence  the  trajectory  record  would  end. 

7 . 5  Particle  instantaneous  velocities 

A  primary  objective  of  the  study  of  axonal  particle 
movement  has  been  to  understand  the  mechanisms  responsible 
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for  this  motion.  The  relative  magnitude  of  the  inertial  and 
viscous  forces  acting  on  a  particle  can  be  estimated  with 
the  Reynolds  number,  which  gives  the  ratio  between  the 
inertial  and  viscous  forces  acting  on  a  body. 


R  =  L VP 
z ' 

(7.1) 


V  -  velocity  of  the  object 
z  -  dynamic  viscosity  of  the  medium 
p-  density  of  the  surrounding  medium 
L  -  character ist ic  length  of  the  object  (particle  diameter) 

For  an  axonal  particle,  typical  values  used  in  the 
literature  (Cooper  and  Smith, 1974;  Copeland, 1976)  are 
L=2*107  m,  V=2*106  m/sec,  jo=1*103  kg/m3,  z=1  Pa  sec.  These 
values  give  a  Reynolds  number  of  R=4*1010,  which  implies 
that  the  inertial  forces  acting  on  a  particle  are  very  small 
when  compared  with  the  viscous  forces.  Unless  there  is  a 
force  acting  on  a  particle  it  will  remain  at  rest,  and  when 
moving,  will  come  to  rest  almost  immediately  should  the 
motive  force  cease.  As  the  viscous  force  on  an  object  is 
directly  proportional  to  its  velocity,  the  instantaneous 
velocity  of  a  particle  is  hypothesized  to  be  directly 
proportional  to,  and  a  measure  of,  the  forces  acting  on  it. 
Since  examination  of  the  instantaneous  velocities  of  a 
considerable  number  of  particles  has  consistently  shown  that 
particle  velocities  vary  smoothly  with  time,  it  is 
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hypothesized  that  the  net  motive  force  behaves  in  a  similar 
manner.  The  temporal  relationship  may  be  related  to  some 
periodicity  (with  a  period  of  more  than  5  seconds)  in  the 
force  generating  mechanism  which  propels  each  particle. 

The  "noisy"  appearance  of  the  instantaneous  velocity 
estimates  as  obtained  with  the  incremental  distance  method 
and  reported  in  the  literature  (Breuer  et  al ., 1975; Forman  et 
al ., 1977a; Lees tma  et  al.,1977)  agrees  closely  with  that 
predicted  by  simulation  analysis  in  this  work.  As  the  errors 
in  instantaneous  velocity  estimates  obtained  with  this 
method  have  been  shown  to  be  considerable,  extensive 
comparison  of  the  instantaneous  velocity  results  and 
associated  parameters  obtained  in  this  work  with  those  found 
in  the  literature  is  not  appropriate. 

Examination  of  the  power  spectra  and  instantaneous 
velocities  both  led  to  the  conclusion  that  the  instantaneous 
velocities  of  somatopetal ly  moving  particles  tended  to  vary 
more  widely  about  their  mean  values  than  did  that  of 
somatofugally  moving  particles.  Results  also  indicated  that 
the  variation  of  a  particle's  instantaneous  velocity  was  not 
related  to  its  mean  value. 

7.6  Spatial  analysis 

Previous  attempts  to  identify  spatial  patterns  in 
particle  motion  have  been  crude.  In  their  analysis  of 
particle  motion  in  neurites  in  tissue  culture,  Breuer' s 
group  plotted  estimates  of  particles'  instantaneous  velocity 
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against  distance  along  the  neurite  in  an  attempt  to  identify 
longitudinally  separated  zones  of  high  or  low  velocity 
(Breuer  et  al.,1975).  Their  results  revealed  no  such 
regions,  and  they  concluded  that  if  zones  of  high  or  low 
velocity  do  exist,  they  do  not  extend  uniformly  across  the 
neurite.  Leestma' s  group,  using  preparations  and  methods 
similar  to  Breuer,  also  concluded  that  zones  of  uniform 
velocity  do  not  extend  across  the  neurite  (Leestma  et 
a  1., 1977).  They  did,  however,  note  different  areas  in  which 
different  particles  exhibited  similar  "disturbances  of 
motion"  at  different  instants  in  time.  This  evidence  led 
them  to  suggest  that  areas  of  disturbance  are  themselves 
moving  along  the  axon. 

It  is  the  opinion  of  this  author  that  these  conclusions 
are  largely  unsubstantiated  when  the  preparations  and 
methods  used  are  considered.  The  neurites  in  tissue  culture 
used  by  both  groups  were  so  small  that  analysis  of  the 
details  of  motion  must  have  been  very  difficult.  Questions 
have  also  been  raised  as  to  whether  behavior  observed  in 
cultured  neurites  can  be  considered  typical  of  that  in 
mature  axons  (Forman  et  al., 1977a).  The  considerable  error 
in  the  instantaneous  velocity  estimation  procedure  used  by 
both  groups  has  been  previously  demonstrated.  Leestma  notes 
when  discussing  their  preparation  that  scanning  electron 
micrographs  of  the  nerve  fibres  filmed  indicated  that  rather 
than  being  single  neurites  they  could  in  fact  be  composed  of 
a  number  of  neurites  bound  closely  together.  Hence  particle 
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trajectories  and  "disturbances  of  motion"  appearing  on  film 
to  be  in  the  same  neurite  could  actually  be  in  a  number  of 
adjacent  neurites.  Assuming  regions  of  uniform  velocity  do 
exist,  it  seems  unlikely  that  correspondingly  located 
similar  regions  will  be  found  in  adjacent  neurites. 

In  the  preparations  used  in  this  work,  all  the 
trajectories  filmed  from  a  given  preparation  occurred  within 
a  single  axon.  In  addition,  the  axons  filmed  were 
considerably  larger  than  the  cultured  neurites  used  by 
Breuer  and  Leestma,  and  as  a  result  more  trajectories  were 
tracked  from  a  given  segment,  with  a  corresponding  increase 
in  the  velocity  estimates  available  for  analysis. 

Overlapping  the  trajectory  plots  of  42  somatopetal  and 
4  somatof uga 1  particles  recorded  from  a  single  segment  of  an 
axon  (figure  4.2)  failed  to  reveal  the  presence  of  tracks 
along  which  the  particles  moved. 

The  spatial  relationship  between  microtubules  and  rod 
shaped  mitochondria  has  been  previously  noted  and  discussed 
in  the  literature,  as  has  the  hypothesis  that  microtubules 
act  as  "tracks"  along  which  particles  move  (Schmi tt , 1969; 
Smith, 1973;  Hammond  and  Smith, 1976;  Smith  and  McLeod, 1979; 
Smith  and  McLeod , unpub 1 i shed  data).  It  has  been  also  shown 
that  a  minimum  level  of  ATP  is  required  to  maintain  fast 
axoplasmic  transport  (Ochs, 1974).  As  mitochondria  are  a 
source  of  intracellular  ATP,  it  has  been  hypothesized  that 
this  energy  source  is  used  to  "power"  particle  motion 
(Schmi tt , 1969) .  To  investigate  the  possible  positional 
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relationship  between  high  particle  velocity  and  mitochondria 
a  transparency  of  the  axonal  structures  and  particle 
trajectories  (figure  4.2)  was  overlaid  on  the  contour  plot 
of  the  velocity-position  data  (figure  6.14).  Visual 
inspection  of  the  two  overlaid  plots  showed  no  apparent 
relationship  between  the  rod-shaped  mitochondria  and  regions 
of  high  particle  velocity. 

The  contour  plots  of  particle  instantaneous  velocity  as 
a  function  of  position  have  indicated  a  relationship  between 
instantaneous  velocities  of  particles  and  their  spatial 
locations.  The  spatial  relationship  may  be  related  to 
factors  such  as  local  variations  in  axoplasmic  viscosity. 

7.7  Conclusions 

This  thesis  has  outlined  the  development  of  a  new 
method  for  the  estimation  of  the  instantaneous  velocity  of 
axonal  particles.  The  superiority  of  the  new  method  over 
that  commonly  used  in  the  literature  was  demonstrated  by 
comparison  of  the  errors  inherent  in  the  use  of  each.  The 
instantaneous  velocities  as  estimated  with  the  new  method  of 
one  hundred  and  seventy-two  particles  were  grouped, 
quantified,  and  compared  using  a  number  of  techniques,  in  an 
attempt  to  identify  patterns  or  relationships  in  particle 
motion.  Although  the  results  were  at  times  inconclusive, 
further  analysis  of  the  instantaneous  velocities  was  beyond 
the  scope  of  this  work. 

The  results  support  the  hypothesis  that  the  movement  of 
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particles  in  axons  is  both  time  and  position  dependent. 
Although  the  mean  velocities  of  the  somatopetal  and 
somatofugal  populations  were  not  significantly  different, 
differences  in  the  average  power  spectra  of  the  two 
populations,  and  in  the  frequency  distributions  of  the 
populations'  instantaneous  velocities  indicate  that  there 
may  be  significant  differences  in  the  motive  mechanisms 
acting  on  each. 

Clearly  considerable  work  remains  to  be  done  in  the 
study  of  axonal  particle  motion.  As  the  work  of  each 
researcher  reveals  additional  details  of  the  underlying 
mechanisms  the  overall  process  may  be  better  understood.  It 
is  the  hope  of  the  author  that  the  methods  developed  in  this 
thesis  will  aid  in  this  endeavor. 
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